Assessing bird habitat occupancy from gradient-based landscape metrics and principal polar spectral indices in the colombian andean region

The patch-mosaic model (PMM) is the most common way to describe the landscape in ecological research. Despite this, the gradient model (GM) was proposed as a more accurate representation of the heterogeneity of landscapes; however, little has been explored on the behavior and performance of continuous variables and surface-based metrics from GM under different analytical scenarios. We address the question: which landscape metrics, patch-based or surface-based, best explain habitat occupancy patterns of six bird species with different ecological preferences? We generated detection histories for six bird species in a fragmented Andean landscape from Colombia. We obtain patch-based metrics from a land cover classification and surface-based metrics from the principal polar spectral indices (PPSi) to describe the landscape. Finally, we fitted dynamic occupancy models using variables derived from landscape models and compared their performance using quasi-AIC for each species. We obtained 909 detections for the six selected bird species. We found that PPSi and surface-based metrics were more informative when assessing occupancy patterns for five of the six species studied. In addition, surface-based metrics allowed to detect interspecific differences between species beyond an affinity for a particular cover type. Surface-based metrics can be an alternative for assessing species response to landscape heterogeneity, particularly those that may be more sensitive to fine-scale changes in vegetation cover. However, there is no single “best” model to describe the landscape for all cases. PPSi can be very useful for land cover analysis in landscape ecology studies as an alternative to more popular vegetation indices.


Introduction
Assessing how species, populations, or communities are related to landscape structure is a central issue in landscape ecology and partly depends on the conceptual model chosen to describe landscape composition and configuration (Lindenmayer et al. 2007;Lausch et al. 2015;Frazier & Kedron 2017). There are different conceptual models to describe landscape structure, although the Patch-Mosaic model (PMM) (Forman 1995) is currently the most popular one (Lausch et al. 2015;Fardila et al. 2017 ). The PMM is an expansion of the patch model, which originates from island biogeography theory (Shafer 1990), and assumes that landscapes are a mosaic of discretely delineated habitat units or homogeneous areas without internal variation (Forman 1995). Despite its popularity and acceptance, the PMM model has been criticized for ignoring or simplifying the continuous nature of environmental gradients (McGarigal et al. 2009;Frazier & Kedron 2017 ). As an alternative, the Gradient Model (GM) has been proposed as a more adequate way of conceptualizing and analyzing landscape structure capable of capturing more heterogeneity and overcoming the limitations of the PMM (McGarigal and Cushman 2005;McGarigal et al. 2009). Thus, an essential task for landscape ecologists is to elucidate the best approach or model to describe the landscape structure on which their research is focused.
In the PMM, the surface describing the landscape typically corresponds to a land cover or vegetation classification (Forman 1995). Under this model, most metrics available to characterize the landscape structure focus on analyzing the spatial pattern of the mosaic of habitat patches or landscape elements (Kupfer, 2012;Lausch et al. 2015;Li & Yang 2015) at three levels: patch, class (several patches of the same type), and landscape (K. McGarigal 2002;Lausch et al. 2015). There are many metrics associated with patch and class level to analyze the spatial pattern of landscape elements. Most of these metrics have a high degree of correlation and quantify patch area and shape, the number or density of edges, degree of isolation, and structure of the surrounding matrix (McGarigal and Marks 1995;McGarigal 2002;Hesselbarth et al. 2019). Despite this, many of these metrics do not represent specific aspects of habitat selection by organisms, and they depend on the spatial scale of analysis and the habitat definition employed (Lausch et al. 2015). Studies of the effects of landscape structure or dynamics on biodiversity typically use metrics based on the PMM (Fardila et al. 2017), mainly because they are easy to obtain and interpret, and require less computational power (Lausch et al. 2015).
In contrast, in the GM, the landscape is a continuous tridimensional surface (x, y, z) where the only discrete unit is the pixel or cell in a raster matrix (McGarigal and Cushman 2005). The x and y raster coordinates indicate the spatial or horizontal location of each pixel, while the pixel value indicates the height (z). If all pixels in the raster have the same value, a homogeneous or flat surface is formed. Studies analyzing landscape structure under the GM have used raster surfaces describing the amount of plant biomass or photosynthetic vigor, such as normalized difference vegetation index (NDVI) (McGarigal et al. 2009) or percent of vegetation cover (Bruton et al. 2015). Unlike the PMM, there are no established guidelines to represent habitat variation in the landscape because the GM has been less explored (McGarigal and Cushman 2005;McGarigal et al. 2009;Lausch et al. 2015). GM landscape metrics, or surface-based metrics, can be categorized into three major groups according to McGarigal et al. (2009): amplitude metrics, surface bearing metrics, and spatial metrics. Amplitude metrics quantify the distribution of height values (i.e., vertical distribution) of the surface (e.g., variability, kurtosis, and skewness). Surface bearing metrics quantify surface texture based on the Abbott-Firestone curve obtained from the inversion of the cumulative histogram of surface height distribution. Spatial metrics measure the horizontal distribution (spatial pattern) of the surface height values.
Analyzing landscape environmental variation using surface metrics has several advantages over patch-based metrics ( McGarigal and Cushman 2005;McGarigal et al. 2009). For example, surface metrics do not require delimiting homogeneous areas or defining habitat edges a priori. Instead, each locality within the landscape is considered a combination of environmental dimensions to which organisms respond, according to their functional traits, ecological preferences, or niche breadth (McGarigal and Cushman 2005;McGarigal et al. 2009). This makes it possible to include all locations within the landscape for analysis, independent of whether we consider them habitat or surrounding matrix, as well as facilitating species-centered modeling and interpretation (i.e., exploring habitat preferences for each species and analyzing responses to fragmentation of that habitat) (Lausch et al. 2015). On the other hand, surface metrics imply managing, processing, and storing large amounts of multidimensional spatially explicit information, that are less intuitive, and still need to be standardized in their use and interpretation (Lausch et al. 2015). Several of these drawbacks are becoming less problematic. There is now a large availability of data to analyze landscapes in different environmental dimensions and spatial scales resulting from remote sensing of the Earth's surface (e.g., Landsat, Sentinel) (Wang & Gamon, 2019). Moreover, using available open-source tools for obtaining, visualizing, processing, and analyzing geographic and remote sensing data has also increased in ecology and conservation sciences (e.g., QGIS, Google Earth Engine, geodiv package for R) (Rocchini et al. 2017;Kumar & Mutanga 2018;Smith et al. 2021).
In this study, we compare the utility of PMM and GM as approaches to describe how bird species perceive habitats at the landscape scale in a fragmented Andean landscape from Colombia. Specifically, we address the question: Which landscape metrics, patch-based or surface-based, best explain habitat occupancy patterns of six bird species with different ecological preferences? We expect surface metrics generated under the GM to be more informative than traditional binary or proportion variables based on the PMM regardless of species ecologies. Nonetheless, for habitat specialists, variables based on PMM may be as informative as those generated under the GM. In addition, considering that there is no standardized continuous surface representing the landscape analogous to vegetation cover classification, we tested the effectiveness of principal polar spectral indices (PPSi) (Moffiet et al., 2010) as an alternative to more traditional metrics (e.g., NDVI) to represent variation in vegetation cover under the GM. The vegetation cover type and its variation across the landscape have been related to aspects of bird ecology such as feeding, dispersal ability, behavior, reproduction, and predation probability (Bélisle et al. 2001;Barlow et al. 2007;Kennedy et al. 2010;Neuschulz et al. 2013;Carrara et al. 2015;Walter et al., 2017;Geoffroy et al., 2019). Thus, we expect that species responses to PPSi indices and greenness-based surface metrics can be interpreted in terms of habitat use preference.

Study region
Our study region was located on the eastern side of the Central Cordillera of Colombia, corresponding to a highly fragmented landscape of the humid tropical biome between 150 and 2500 m elevation (Fig. 1). This region covers a mixture of sub-Andean humid forest relicts of different sizes, agroecosystems, shrublands, secondary vegetation in different stages of succession, cattle pastures, and artificialized soils (IDEAM et al., 2017). Also, this region contains three large artificial water bodies corresponding to the San Lorenzo, Playas, and Punchiná hydroelectric dams.

Bird sampling
We chose six bird species commonly sighted in the study region with different habitat requirements and land-use patterns: Formicarius analis, Oncostoma olivaceum, Cercomacroides tyrannina, Ortalis columbiana, Troglodytes aedon and Tyrannus melancholicus. We generated detection histories for each selected bird species by establishing 70 fixed sampling points along 14 transects (1.25 km each) covering the transition from forest to secondary vegetation or grassland (Fig. 1). A minimum of 215 linear meters separated each observation point. We visited each sampling point during June and July of 2014 and 2015, recording all visual or auditory detections for 15 min within a maximum radius of 100 m. Bird records were taken by at least four experienced observers, visiting each site for two consecutive days to obtain four temporal replicates per year. All sampling events were conducted between 6:00 am, and 10:00 am.
F. analis (Formicariidae) is a medium ground species (body mass: 62.19 g; length: 19 cm) that feeds mainly on invertebrates found on the forest floor and is distributed in Central and South America. O. olivaceum (Tyrannidae) is a small insectivorous species (body mass: 6.6 g; length: 9.1 cm) found in forest understory, secondary vegetation, and shrublands areas; this species is distributed in Panama and northern Colombia. C. tyrannina (Thamnophilidae) is a small insectivorous species (body mass: 9.34 g; length: 16 cm) that feeds mainly on invertebrates, and is commonly observed in forests, forest edges, and secondary vegetation. O. columbiana is a large species endemic to the Andean region of Colombia (body mass: 600 g; length: 53 cm) that feeds mainly on fruits and seeds. This species is common in its distribution range and is associated with forest edges, gallery forests, forest clearings, and pastures with scattered trees; it is occasionally observed in urban and peri-urban areas. T. aedon is a small species (body mass: 10.85 g; length: 11.4 cm) that inhabits semi-open areas, forest clearings, peri-urban and urban areas, and feeds on invertebrates and fruits. Finally, T. melancholicus is one of the most widely distributed species in Colombia, inhabiting pastures, grasslands with scattered trees, urban and peri-urban areas. It is a medium sized species (body mass: 35 g; length: 22 cm) with generalist feeding habits

Land cover classification and Patch-based metrics
Our study required a detailed digital land cover map, which we generated from 12 scenes of the MSI (Multispectral Imager) sensor in the Sentinel 2 Earth observation satellite, employing the Sen2Cor and Sen2Three modules of the free Sentinel Application Platform (SNAP) software (Louis et al. 2016;Main-Knorn et al., 2017). The Sentinel 2 MSI imagery spatial resolution varies from 10 to 60 m per pixel, depending on the spectral band (European Space Agency -ESA, 2015); in this case, the spatial resolution of the performed analyses was 10 m. The combination algorithm gives greater importance to suitable pixels (i.e., no clouds or shadows) from scenes temporally closer to 2015, using pixels from more distant scenes only if closer scenes had low-quality information. Using the Sentinel 2 imagery, the Continental, coastal, and marine ecosystems map of Colombia (IDEAM et al., 2017), and a high-resolution aerial photograph, we performed a supervised land cover classification using the open-source Semi-Automatic Classification Plugin (SCP) for QGIS 3 (Congedo 2016), identifying five vegetation -land cover types: forest, secondary vegetation, grassland, water and bare or urban soil (Fig. 2). The obtained Kappa index was 0.82, and the overall supervised classification accuracy across the landscape was 88.07%, indicating an acceptable performance (Olofsson et al., 2014) (see supplementary material Table S2).
We calculated 24 patch-based metrics for each sampling unit focusing on forest cover using the landscapemetrics package for R (Hesselbarth et al. 2019). Since many of these metrics are strongly correlated with each other, we retained only those metrics that showed the lowest degree of correlation (i.e., Pearson r < 0.7) with any of the other metrics: total area (Ca), average shape index (Shape), number of patches (Np), clumpiness index (Clumpy), and mean Euclidean distance to nearest neighbor (Enn). The Ca and Shape metrics quantify patch area and shape, while Np, Clumpy, and Enn are aggregation metrics. We implemented circular nested buffers with four different radii (i.e., 70, 130, 270, and 8100 m) to calculate the metrics and obtain a multiscale average for each sampling point. The equations and descriptions of each patch-based metric used here are widely described in the scientific literature (e.g., (McGarigal and Marks 1995).
Generally, Ca expresses the average area, in hectares, of all forest patches within the assessed radii and is equivalent to the habitat amount according to (Fahrig 2013), assuming forest cover as a proxy for habitat of forest specialist species. The average Shape index expresses the ratio between perimeter and area of the forest patches located within the assessed radii. The shape index shows low values for regular-shaped patches (i.e., quadrangular). The number of patches, Np, is one of the simplest ways to quantify fragmentation in a landscape; highly fragmented landscapes have many patches. Clumpy indicates the aggregation degree of forest patches and ranges from − 1: maximum disaggregation to 1: maximum aggregation. Finally, Enn measures the average edge-to-edge distance to the nearest neighbor of the same class, in this case, forest. We expect Ca, Clumpy and Enn to be important for habitat specialists such as Formicarius analis, Oncostoma olivaceum and Cercomacroides tyrannina, whereas the number of patches and clumpy to be informative for generalist species such as Tyrannus melancholicus and Ortalis columbiana.

Land cover gradient representation and surface-based metrics
We obtained the principal polar spectral indices (PPSi) from the polar transformation of the Sentinel 2 spectral bands for the study region, adapting the methodology developed by (Moffiet et al., 2010) for ETM + and TM sensor scenes from the Landsat mission (See Supplementary material Fig. S1). The PPSi indices decompose the land cover variation into three orthogonal dimensions: greenness (PPSG), brightness (PPSB), and wetness (PPSW), similar to the Tasseled Cap indices (Kauth and Thomas 1976). With PPSi indices, it is possible to differentiate the variation in land cover between sites associated with changes in the vegetation itself (i.e., type of cover, amount of foliage, proportion of vegetation) from the variation associated with spectral brightness or the presence of water, unlike unidimensional indices such as NDVI (Moffiet et al., 2010). The PPSW index identifies the presence of water. The PPSG index represents the most significant variation among the vegetation cover types classified for the study region, separating forest and secondary vegetation from pasture and bare or urbanized soil areas. Therefore, using the PPSB index in conjunction with the PPSG index is necessary to separate forest and secondary vegetation cover. Thus, sites or pixels with high PPSG values will have high PPSB values if they correspond to secondary vegetation and low PPSB values if they correspond to a forest (See Supplementary material Fig. S1). Combining these three indices allows describing land cover through continuous variables that accurately represent the vegetation gradient present in a region based on information derived from remote sensing (Fig. 2).
We calculated 11 surface metrics for each sampling unit based on the PPSG greenness index using the geodiv package for R (Smith et al. 2021). We calculated surface metrics by implementing circular buffers with the same radii used to calculate patch-based metrics (i.e., 70, 130, 270, and 8100 m) to obtain a multiscale average per sampling point. Then, we retained five metrics that showed a lower degree of correlation with each other (i.e., Pearson r < 0.7): average roughness (Sa), surface asymmetry (Ssk), surface kurtosis (Sku), surface area ratio (Sdr), and summit density (Sds). Sa, Ssk and Sku are amplitude metrics quantifying the variability and statistical distribution in the surface values (i.e., vertical distribution) without considering the spatial arrangement, location, or distribution of the peaks and valleys formed by the surface values (i.e., horizontal distribution). Sdr and Sds are spatial metrics, i.e., they consider both vertical and horizontal distribution of surface values. The interpretation of surface-based metrics in landscape ecology depends on the variable or habitat feature that is represented on the raster surface.
In this study, Sa measures the absolute deviation of the PPSG values from the mean within each radius and can be interpreted as the vegetation cover variation; high Sa values indicate high heterogeneity in vegetation cover. Sa increased in sites composed of contrasting vegetation cover in our study area, such as forest-grassland edge zones with forest fragments of various sizes surrounded by different cover types. Ssk measures the skewness, and Sku measures the kurtosis of the distribution of PPSG values, so they are interpreted as complementary measures of dominance. Low Ssk values indicate dominance of high greenness values (i.e., forest or secondary vegetation), and high Sku values are interpreted as high dominance of greenness values around the mean, i.e., there is a low probability of obtaining values far from the mean. Sku is not a directional metric, so it does not indicate the dominant greenness value specifically in each locality. Sku may increase in relatively homogeneous sites around the average greenness value, which may be low or high. In our study area, Sku correlates moderately with greenness (PPSG; r = 0.48, p < 0.001), summit density (Sds; r = 0.48, p < 0.001) and forest area cover (Ca; r = 0.46, p < 0.001), indicating that higher Sku values tends to occur in forest or secondary vegetation areas predominantly.
The area surface ratio, Sdr, measures the ratio between the observed surface area and the area of a flat surface with the exact x,y dimensions, so it increases with increasing local variability in slopes; Sdr is an indicator of the amount and magnitude of vegetation cover contrasts or edges. Sds measures the number of local peaks of greenness per area and is therefore related to the number of pixels/sites with high greenness values (forests and secondary vegetation) and their spatial distribution. Sds can be interpreted as an indicator of forest fragmentation since the density of greenness peaks (i.e., localities with high PPSG values) corresponds to forest or secondary vegetation cover patches. Highly fragmented landscapes will have low Sds values, i.e., low pixel density with high PPSG values.

Species occupancy analysis
We fit dynamic occupancy models (MacKenzie et al. 2003;Kéry & Royle 2021) for the six selected bird species using metrics from the PPM and GM landscape models as occupancy covariates (Ψ) and the date and sampling time as detection covariates (p). No covariates were used to predict colonization (γ) or extinction (ε), assuming little or no variation in the dynamics governing the individual's distribution of the species of interest over this period (i.e., one year) (Betancur et al. 2020). We used the colext function in the unmarked package for R software to fit all models (Fiske and Chandler 2011). First, we built a global model for each species with all selected landscape metrics plus elevation to test the general model fit and diagnose the presence of overdispersion with the mb.gof.test function from the AICcmodavg package. We tested the goodness-of-fit of the global models by comparing the observed chi-square statistic with the reference distribution generated by 999 simulations. We evaluated overdispersion by calculating the variance inflation factor (cˆ). Global model is theoretically the best-fitting model because all candidate models are special cases of the global model, so we can compute from it an unambiguous cˆ to diagnose the presence of overdispersion (Burnhan and Anderson 2002).
We modeled the detection parameter for each species while holding occupancy, colonization, and extinction constant. We then modeled the occupancy parameter using the top-ranked detection sub-model and holding colonization and extinction constant. We built a set of 42 candidate models for the occupancy parameter avoiding models with correlated covariates or with more than three covariates to avoid problems of model convergence and interpretation of results (See Supplementary material Table S3). Our main objective was to compare the metrics derived from the PMM and GM landscape models, so we did not include models combining patch-based and surfacebased variables, except for the global model. Models were evaluated according to the quasi-AIC criterion (QAIC), calculated using the global model variance inflation factor. Selecting the best models from the QAIC and cˆ allows considering the overdispersion in the calculation of the model parameter standard errors caused by possible lack of independence in the predictors or species detection histories (Richards, 2008). We ranked models using ΔQAICc and considered that models with ΔQAICc < 2 had the strongest support (See Supplementary material Table S4). For the final selection, we excluded models with poor explanatory power or uninformative variables, i.e., variables whose confidence interval included zero (Arnold, 2010). Finally, we assessed the goodness-offit of the final models for each species using the mb. gof.test function on the AICcmodavg package. With the final models for each species, we made inferences about the covariates' effects and the direction of the relationships.

Results
We obtained 909 detections for the six selected bird species based on 25.5 effective hours of sampling (i.e., 15 min x 70 sites x 8 replicates). Ortalis columbiana had the highest number of detections and naïve occupancy (Ψ naïve = 0.66), while Troglodytes aedon had the lowest number of detections (Table 1). The detections varied from 1.5% (Tyrannus melancholicus) to 43% (Cercomacroides tyrannina) of all detections between the two sampling years. We obtained a single final top-ranked dynamic occupancy model for each bird species after our selection procedure and removal of models with uninformative covariables. According to goodness-of-fit tests, the final top-ranked models for each species performed well (Table 2). According to the top-ranked models, the occupancy probabilities of the assessed bird species were primarily associated with continuous variables and surface-based metrics, except for O. columbiana (Table 2).
Occupancy models revealed that bird species respond to different landscape features related to vegetation cover variability, distribution, and isolation  Fig. 3). The occupancy probability of F. analis was positively associated with summit density (Sds), indicating a habitat use preference for localities dominated by high greenness values, i.e., continuous, or little fragmented forest vegetation. According to the top-ranked model, the occupancy probability of F. analis increased (Ψ > 0.75) in localities surrounded by the greatest amount of forest in our study area (Sds > 0.05). The occupancy probability of C. tyrannina increased in sites with low surface-area ratio values (Sdr), indicating a preference for areas with no abrupt changes in greenness, i.e., sites with continuous cover and away from forest edges.
Finally, O. columbiana was associated with a patch-based metric measuring the degree of forest cover fragmentation or subdivision. The occupancy probability of O. columbiana was best explained by elevation and the number of forest patches (np). The model obtained indicated that this species prefers sites above 1000 m and highly fragmented forest cover.

Discussion
The gradient model (GM) was proposed as a more accurate representation of the natural heterogeneity of landscapes; however, little has been explored on the behavior and performance of continuous variables and surface metrics under different analytical scenarios (McGarigal et al. 2009;Kedron et al. 2018). Our work advanced this topic by using metrics related to GM to analyze occupancy patterns of six bird species with different ecological requirements compared to traditional patch model metrics (PMM). We found that GM-based metrics were informative when assessing occupancy patterns in response to landscape land cover structure for most of the species studied (i.e., the top-ranked occupancy models included GM covariates instead PMM covariates). The increasing availability of spatially explicit data (i.e., remote sensing) opens the possibilities of using these variables and the GM in many landscape applications.
We expected patch-related variables to be informative for bird species with forest affinity. On the other hand, all forest-dependent species responded mainly to surface-based metrics, i.e., the occupancy of these species in each location depends on both the vegetation cover of the site and the surrounding cover in the local landscape. In general, these species responded negatively to increased variability in cover (i.e., surface roughness [sa]), preferring areas dominated by high greenness values (i.e., high values of kurtosis [sku] and skewness [skk], and low values of surfacearea ratio [sdr]). This indicates that forest species are sensitive to subtle changes in vegetation cover, which are not quantified under the patch model, and may be related to low tolerances to habitat conditions modified by habitat degradation or edge effects, i.e., lower niche breadth relative to the dimensions represented in vegetation covers. This is characteristic of the forests found in the study area, which have been disturbed in a variety of ways, resulting in a heterogeneous mixture of secondary forests that are currently under protection. Similar results have been found in other studies highlighting the importance of local vegetation characteristics and not just the landscape configuration or composition (Morante-Filho et al. 2021). The relationships found may also be related to a low capacity for movement or dispersion, impeding the search for and attaining resources in variable environments. For example, F. analis -a terrestrial antbird -and C. tyrannina -an understory antbird -that both occasionally follow army ants, may require good quality vegetation with a rich understory to move (Hilty and Brown 2001).
We expected surface-based metrics to be more informative for generalist bird species associated with a variety of habitats. This was the case of Tyrannus melancholicus where both polar indices were informative, and Troglodytes aedon, but not the case for Ortalis columbiana. T. melancholicus was found to associate with open habitats, nonetheless, other studies have found that this can be an artifact of detection bias (Ruiz-Gutiérrez et al. 2010). T. aedon, a species associated to human settlements did not exhibit any preferences for a specific cover type, which is consistent with the idea that most of the landscape is perturbed and suitable for the persistence of this species. Only the Colombian chachalaca, O. columbiana, was positively related to the number of forest patches (np), a PMM covariate, indicating a preference for highly fragmented landscapes. In our analysis, none of the GM-based metrics captured the pattern of forest fragmentation in the same way that the np metric did. The Colombian chachalacas often use forest patches or edges to perch and vocalize in groups, moving between patches to feed or to escape. Although, this species is not a forest specialist, the spatial configuration of forest patches may be important to their occupancy patterns and habitat use (Borges 1999;Rivera-Gutiérrez 2006). We calculated GM-based metrics using the PPSG index, which by itself does not allow differentiation between forest and secondary vegetation cover, which may have influenced the result obtained for this species.
Considering our results, we adopt the view of Lindenmayer et al. 2007 supporting a pluralistic view that recognizes the importance of both continuous and categorical variables in the analysis of landscape variation. Studies explicitly comparing the performance of both PMM and GM models have generally concluded that it is necessary to adopt a pluralistic view in landscape analysis (Lindenmayer et al. 2007), deriving variables or metrics from different landscape models, as the performance of these depends on the behavior or habitat use of species and the characteristics of the landscape under study (Price et al. 2009;Bruton et al. 2015;Salgueiro et al. 2018). For example, (Salgueiro et al. 2018) found that the GM performs better in describing the landscape for specialist bird species, which may be more sensitive to slight variations in landscape characteristics. (Price et al. 2009) found that savanna bird species in Australia that take advantage of different vegetation cover types for foraging or nesting fit the GM better, while species that nest in dense vegetation areas exclusively fit the PMM model better.
Although at a broad level the assessed species could be classified as forest-specialists or generalists, with a preference for areas with little woody vegetation, the GM allowed us to detect interspecific differences for the understanding of their ecology and conservation. For example, forest species such as O. olivaceum, C. tyrannina and F. analis were related to different variability metrics and spatial distribution of greenness values. O. olivaceum occupancy was positively related to skewness (Sku), meaning preference for sites (or local landscapes) highly dominated by forest and secondary vegetation. C. tyrannina was inversely related to surface -area ratio (Sdr), indicating preference for areas with continuous forest cover away from edges, and F. analis was positively associated with summit density (Sds), meaning preference for sites with a low fragmentation degree. Detecting these types of differences can mean progress in proposing specific conservation measures for a particular species or understanding co-occurrence patterns between species beyond an affinity for a particular cover type (Abdel Moniem and Holland 2013). Additionally, including site values of greenness (PPSG) or brightness (PPSB) and surface-based metrics in analytical models can help us to estimate whether species respond to the presence of a cover type at a site and whether they depend on the existing cover at surrounding sites (i.e., local landscape).
To ensure an adequate description of the landscape with both PMM and GM, we must obtain a surface representing habitat composition (e.g., vegetation covers, NDVI, % vegetation cover, PPSi) and then estimate landscape metrics describing their vertical and horizontal variation (e.g., patch-based, or surfacebased metrics) (McGarigal and Cushman 2005; Frazier & Kedron 2017). Selecting surface and landscape metrics should be based on knowledge of the habitat requirements and behavior of the species assessed (Schindler et al. 2015); however, we rarely have relevant information on these aspects of species ecology, especially in the tropics. Because of this, a cost-efficient alternative is to use vegetation cover as an indicator for habitat and obtain diverse metrics to quantify variability pattern expecting that some of these metrics or a combination of them will be informative for most species (Frazier and Kedron 2017). Considering that in our analysis the landscape metrics derived under either landscape model come from the same base information (i.e., Sentinel 2 imagery) and represent the same landscape characteristics (i.e., vegetation cover), the difference in the occupancy models performance and the conclusions that can be drawn from them are interesting. This means that the limitation of the patch-based landscape metrics for most of the species analyzed may not necessarily be due to vegetation cover being a poor indicator of the habitat characteristics to which species respond but to the loss of information associated with its discretization, homogenization, and simplification (McGarigal and Cushman 2005;McGarigal et al. 2009). In this sense, using principal polar spectral indices (PPSi) as an indicator of vegetation or land cover in combination with surface metrics can improve our understanding of the biodiversity patterns at the landscape scale, especially in highly heterogeneous landscapes such as the ones included in this study and highly frequent for tropical latitudes.
The principal polar spectral indices (PPSi) have two significant advantages over other and more popular forms of describing vegetation cover in a continuous format, such as percent tree cover or normalized difference vegetation index (NDVI) (Moffiet et al. 2010). First, from an information standpoint, PPSi indices are more sensitive to slight differences in vegetation cover composition among sites, even those classified with equal cover in the patch model (Moffiet et al. 2010;Ramdani et al. 2019). Second, the greenness index (PPSG) directly quantifies aspects associated with vegetation and allows differentiating major cover types, making it possible to interpret surface metrics in terms of variability, dominance, or contrast in vegetation cover types (Moffiet et al. 2010). As cover type is a popular indicator of habitats in landscape ecology field, using surfacebased metrics to describe the vertical and horizontal variation of PPSi indices can help us understand and communicate the results obtained using terminology, ideas, and concepts more familiar to most researchers currently using the PMM and discrete land cover maps. This last point could be crucial considering that one of the most significant limitations to the GM popularization has been how unintuitive surfacebased metrics are compared to patch-based metrics, as many of these originated to evaluate manufactured surfaces with a focus on quality control in the industrial field (Abbot and Firestone 1933;Thomas 1981;McGarigal et al. 2009;Kedron et al. 2018). Exploring behavior of the PPSi indices in different biomes or regions could allow us to test applicability to vegetation mapping in other landscape types (e.g., urban, coastal, savanna, etc.).
According to (Kedron et al. 2018), there are no natural analogs between surface-based and patchbased metrics, so the focus should be directed at conceptualizing surface-based metrics and not using them as indicators of patch-based metrics. One way forward would be to hypothesize associations between surface metrics and vegetation structure based on knowledge of each species natural history, and then evaluate these ideas with on-site measurements. For example, we hypothesize that occupancy of C. tyrannina increases at lower values of Sdr for the PPSG index, indicative of preference for sites with continuous vegetation cover and away from forest edges, and occupancy of F. analis increases at high values of Sds, indicative of continuous forest patches. We can measure these specific vegetation structure properties on site and confirm these associations. One potential complication is the generality of these interpretations for other landscapes.
Finally, our results suggest that using ecologically relevant environmental surfaces for species combined with surface-based metrics that quantify vertical and horizontal variability can generate enough information to characterize ecologically diverse species' responses, particularly those that may be more sensitive to fine-scale changes in vegetation cover. Gradient-based landscape metrics were more informative than patch-based metrics for characterizing the occupancy patterns of most bird species assessed; however, the idea that there is no single "best" model to describe the landscape for all cases or all species persists (Lindenmayer et al. 2007). The GM offers a great opportunity to take advantage of all the high spatial, temporal, and spectral resolution data from remote sensing without losing information due to the simplification or reduction of the categorical classifications of land cover. Therefore, efforts should be made to obtain ecologically informative continuous variables and characterize surface-based metrics' behavior. In this sense, the polar principal spectral indices (PPS) can be very useful as a starting point for land cover characterization, taking advantage of a large amount of information from remote sensing and helping to popularize the GM by allowing its interpretation in the traditional terms and concepts generated from the PPM and the use of categorized vegetation cover maps.