Spatial and Temporal Dependency of NDVI Satellite Imagery in Predicting Bird Diversity over France

: Continuous-based predictors of habitat characteristics derived from satellite imagery are increasingly used in species distribution models (SDM). This is especially the case of Normalized Difference Vegetation Index (NDVI) which provides estimates of vegetation productivity and heterogeneity. However, when NDVI predictors are incorporated into SDM, synchrony between biological observations and image acquisition must be questionned. Due to seasonal variations of NDVI during the year, landscape patterns of habitats are revealed differently from one date to another leading to variations in models’ performance. In this paper, we investigated the inﬂuence of acquisition time period of NDVI to explain and predict bird community patterns over France. We examined if the NDVI acquisition period that best ﬁt the bird data depends on the dominant land cover context. We also compared models based on single time period of NDVI with one model built from the Dynamic Habitat Index (DHI) components which summarize variations in vegetation phenology throughout the year from the fraction of radiation absorbed by the canopy (fPAR). Bird species richness was calculated as response variable for 759 plots of 4 km 2 from the French Breeding Bird Survey. Bird specialists and generalists to habitat were considered. NDVI and DHI predictors were both derived from MODIS products. For NDVI, ﬁve time periods in 2010 were compared, from late winter to begin of autumn. A climate predictor was also used and Generalized Additive Models were ﬁtted to explain and predict bird species richness. Results showed that NDVI-based proxies of dominant habitat identity and spatial heterogeneity explain more bird community patterns than DHI-based proxies of annual productivity and seasonnality. We also found that models’ performance was both time and context-dependent, varying according to the bird groups. In general, best time period of NDVI did not match with the acquisition period of bird data because in case of synchrony, differences in habitats are less pronounced. These ﬁndings suggest that the most powerful approach to estimate bird community patterns is the simplest one. It only requires NDVI predictors from a single appropriate time period, in addition to climate, which makes the approach very operational.


Introduction
Predicting changes in bird species distributions and community patterns is a major challenge in the context of global changes. Over continental and national extents, macro-climatic factors such as temperature and precipitation gradients are the main determinants of species assemblages [1][2][3].
At finer spatial scales, the components of landscape mosaics play a key role in explaining species habitat selection [4][5][6]. Remote sensing-based predictors of habitat characteristics may contribute to improve the performances of species distribution models (SDM) [7]. However, the question of the most appropriate data and representation (discrete vs. continuous-based metrics) in SDM still remains [8][9][10][11][12].
Using unclassified satellite images with continuous-based metrics is a fast way to estimate landscape productivity and heterogeneity. These indirect surrogates can explain and predict community patterns with equal or higher performances than those obtained with land cover maps [11,[13][14][15]. The normalized difference vegetation index (NDVI) which is calculated from near-infrared and red bands is widely used to reflect environmental factors [16]. At large scale, NDVI is a good estimate of ecosystem productivity, enabling to model energy-abundance relationships or to predict species richness [17][18][19]. Locally, its spatial heterogeneity can be related to the diversity of ecological habitats, also explaining patterns of species groups with different ecological preferences [20,21].
Because NDVI varies during the year according to the vegetation phenological stages, this indicator conveys differently the landscape mosaic between the seasons, leading to contrasting performances in explaining biodiversity patterns [21][22][23]. In [21], the authors found that autumn NDVI is the best period to explain and predict species richness of different breeding bird groups in a region of southwestern France. In this area made of heterogeneous landscapes, the unclassified NDVI acquired in autumn allow a better discrimination between harvested crops (low NDVI), woodlands, hedgerows, and grasslands (high NDVI). This discrimination is required to explain bird behaviors of habitat selection or avoidance, according to their nesting and diet preferences. In this study, we aim at extending this result at national scale by assessing how the the NDVI acquisition period that best model the relationships with bird data varies according to landscape types and their dominant land covers.
Instead of using single-date data to explain patterns of biodiversity, recent approaches proposed to explicitely take into account variations in vegetation phenology throughout the year. Among these ones, the Dynamic Habitat Index (DHI) incorporates three components based on the temporal variation of canopy light absorbance through the year including the annual sum, the minimum and the monthly variation [24,25]. This indicator has been already taken up by studies on different taxonomic groups [26][27][28] and a recent study has compared DHI based on different MODIS products to explain breeding bird species richness in North America [29]. However, until now, the added value of this integrative indicator of intra-annual variation, compared to data acquired at a single date, has never been tested.
In this study, we used a large bird dataset collected at the French national scale to (i) compare single date NDVI data and DHI to explain and predict species richness patterns of different bird groups; (ii) identify the best NDVI acquisition date to explain and predict bird species richness; and (iii) test whether the best NDVI acquisition date depends on the dominant land cover in landscapes.

Bird Data
Bird data were obtained from the French Breeding Bird Survey (FBBS), a standardized monitoring program in which volunteer ornithologists identify and record breeding birds year after year [30]. In the FBBS program, each observer provides a point location, and a randomly 2 km × 2 km survey square shaped plot is selected around a radius of 10 km from this location. Each plot contain 10 point counts with 100 m radii which are separated from each other by 300 m at a minimum and scattered in the landscape in order to best represent the diversity of habitats ( Figure 1). Each point count is performed twice a year during the breeding season (from April to mid-June) to record early singing species as well as last migrants. In each point count, heard or seen species are recorded during 5 min between 1 h and 4 h after sunrise. In our study, we focused the analyses on the year 2010 for which the sampling effort was high and well distributed at the national scale (a total of 1091 plots). For each point count in the surveyed plot squares, the maximum count over the two visits was retained for each species. We excluded raptors that have home range sizes (several tens of hectares) much higher than point count areas. We also excluded urban and wetland species because their habitats were poorly sampled compared to agricultural and woodland habitats. Moreover, the DHI components based on fPAR is is not available in the heart of the cities.
Plots were filtered according to the quality of the satellite images available on their locations (see Section 2.2) and 759 plots were finally considered in the analyses. Because birds are expected to respond differently to landscape habitats according to their ecological requirements, four metrics of species richness (SR) per plot square were computed from the species observed at point counts. The first two metrics of SR were related to woodland species and farmland species (specialists of woodland and farmland habitats). The others were related to habitat generalists, and total species richness computed as the sum of all species (see Appendix A for species names). These bird groups (guilds) were based on species information available at the national level (http://vigienature.mnhn.fr).

Climate and Satellite Data
Because avifauna patterns are driven by large climate gradients at the French national scale, a map of macro-climate regions of France was used to derive a climate factor. The delineation of the climate regions in this map is based on 14 temperature and precipitation variables calculated from a 30-year period of climate data  provided by the french national meteorological service [31]. The map is composed of eight regions ( Figure 2): the climate of mountain areas (type 1); the semi-continental climate and the climate of peripheral mountain areas (type 2); the degraded oceanic climate of the north and center plains (type 3); the altered oceanic climate (type 4); the marked oceanic climate (type 5); the altered Mediterranean climate (type 6); the southwest basin climate (type 7); and the marked Mediterranean climate (type 8). Based on this map, the dominant climate region was estimated for each plot of bird data and defined as one qualitative predictor. An alternative would have been to use WorldClim2 data which are locally defined at 1 km 2 scale [32]. However, these data are based on an interpolation using fewer weather stations than the map defined by [31] at French scale. This is why we retained the first source of data (i.e., map of macro-climate regions).
Normalized Difference Vegetation Index (NDVI) was used to provide surrogate of local vegetation cover and landscape heterogeneity patterns over France. The Terra-MODIS NDVI 16-day-composit product at 250-m spatial resolution (MOD13Q1) was retained. To accommodate for seasonal variation and evaluate the temporal effect of NDVI on models performance, five time periods in 2010 were selected ( Figure 2): late winter (6 March 2010, NDVI 065 ), before leaves growing of deciduous tree; late spring (26 June 2010, NDVI 177 ), before winter crops harvesting; late summer (29 August 2010, NDVI 241 ), after winter crops harvesting before summer crops harvesting; end summer (14 September 2010, NDVI 257 ), after partial summer crops harvesting; begin autumn (30 September 2010, NDVI 273 ), after full crops harvesting. At these time periods, landscape composition and heterogeneity are revealed differently in NDVI data. In winter, NDVI values are low in the deciduous forests and in the non-cultivated agricultural parcels because of absence of photosynthetic activity. In spring and summer, NDVI reaches its highest values, mainly in winter and summer crops respectively, but also in grasslands and forests. In autumn, after all crops harvesting, NDVI discriminates well forests and open habitats because of higher difference of NDVI values, agricultural parcels being bare (see Appendix E).
The NDVI data of each time period were selected after checking the MODIS quality layers of all the 16-day time series products of 2010. Based on the average value of the MOD13Q1 pixel reliability over France (0 = good, 1 = marginal data, 2 = snow/ice, 3 = cloudy), the images showing the highest quality (i.e., having the lowest average quality index) were retained for each time period. Then, a filter was applied on the bird survey plots available (n = 1091). Plots including NDVI pixels of bad quality were removed. Only plots composed of 100% of good pixels (i.e., with a pixel reliability = 0) at the five time periods were used (n = 759). Finally, NDVI values in these plots were summarized through several statistical predictors: minimum, maximum, range, average, median, and standard deviation. In order to prevent multicollinearity between covariates, only the average (Avg) and standard deviation (SD) predictors were retained (Spearman's |Rho|< 0.7) [33]. The average predictor is a proxy of landscape productivity and the dominant land cover in the plot. The standard deviation reflects landscape heterogeneity. To examine whether seasonal variations effect of NDVI on the models performance depends on the landscape context, we used land cover/land use datasets to calculate landscape composition in each bird plot. The forest layer was obtained from a nationwide highly detailed vector spatial database (BDTopo ® , French mapping agency IGN) at 1-m resolution. For agricultural areas, the graphic parcel register (RPG) of 2010 was used. The French RPG database is produced annually since 2002 in application of the EU directives related to the Common Agricultural Policy. The land use information, based on a 28 class nomenclature (including several types of crops, as well as artificial and natural grasslands), is obtained from the farmers' declaration of their cultivated areas. The information is defined at the agricultural block level composed of one single field or several contiguous ones. For each block, the proportion of each crop type is provided but without the spatial distribution of crops within the block.
The Dynamic Habitat Index of 2010 was also calculated at the national scale, in order to compare the performance of the DHI-based model with the ones based on single date NDVI. The DHI was derived from the most detailed MODIS fPAR (fraction of Photosynthetically Active Radiation) 8-day composite product at 500-m spatial resolution (MOD15A2H) [34]. It includes three components [26]: (1) the cumulative annual fPAR (DHI cum) reflecting the overall potential of vegetation productivity; (2) the minimum annual fPAR value (DHI min) indicating the minimum amount of vegetated cover; and (3) the coefficient of variation of annual fPAR (DHI cv) expressing the seasonality in productivity ( Figure 3). The average and standard deviation of the three DHI components were calculated for each plot of bird data (n = 759) and used as dependent variables. All the variables were retained in the bird-habitat models after checking correlations between them (Spearman's |Rho|< 0.7). In order to avoid bias in the comparison with NDVI-based models because of difference in spatial resolution, NDVI variables were also computed from the MODIS product at 500-m (MOD13A1). To evaluate the seasonal effect of NDVI on models's performance independently of DHI (see above), we kept the most detailed NDVI data at 250-m. This spatial resolution is likely to produce models with higher explanatory and predictive performances [21].

Statistical Analysis
In a first time, we built models at national scale using the total filtered bird dataset (n = 759). We compared explanatory and predictive performances of NDVI data at several time periods with DHI data. For each of the four species richness metrics, we first constructed a model with only the climate variable, i.e., the dominant climate region of the bird plot, as predictor (e.g., SR woodland~C limate). Then, we added into the climate model both NDVI variables (Avg and SD) of each time period in five separate models (e.g., SR woodland~C limate + Avg.NDVI 065 + SD.NDVI 065 ). Finally, we fitted a final model including DHI variables (6 predictors related to the Avg and SD of the three DHI components, see above) and the climate predictor.
In a second time, we examined whether the best time period of NDVI depended on the plot land cover context. This analysis was only applied to woodland and farmland species richness that were the best explained bird metrics at national scale. We partitioned bird national data into subsets, according to the dominant land cover of the bird plot squares. Landscape composition in each plot was calculated as the proportion of area occupied by four land covers: woody vegetation (including hedgerows), permanent grasslands, winter crops (mainly wheat and barley) and summer crops (mainly maize and sunflower). Because agricultural blocks of the land cover map (RPG) may cover partially the bird plot squares, only the blocks included by at least 70% in the plots were conserved for the estimation. Based on the land cover proportions, four subsets of bird plots were defined with the constraint to obtain a sample size large enough (i.e., with n > 100) and a high dominance of a specific land cover (Table 1). For each subset, we built models with climate and NDVI variables of each date as for the national scale. Table 1. Distribution of the landscape composition in the four subsets of bird plots dominated by woodlands, grasslands, winter crops and summer crops. The average proportion is given with its range (min-max) and standard deviation (SD). Relationships between bird species richness and predictors were explored using generalized additive models (GAMs) with R version 3.1.0 software [35] and the gam package [36]. The general expression of GAM model is [37]:

Plots Dominated by
where Y is the response variable, X i the predictors, f i the non-parametric smooth functions and g the link function relating to the expected value of Y.
GAMs present the advantage to consider non-linear effects between response and predictor variables. Poisson distribution and a log link function were used to calibrate the models. The most parsimonious models (i.e., with the lowest Akaike Information Criterion values) were selected by stepwise backward and forward variable selection [38]. In addition, the procedure was allowed to select from one to four degrees of smoothing (1 = no smoothing) for each quantitative variable in GAMs.
The goodness of fit of the models (explanatory performance) was quantified by the amount of explained deviance (%D 2 ). The independent contribution of each variable to explain species richness was also estimated from a hierarchical partitioning procedure [39].
The predictive performance was calculated by k-fold cross-validation [40]. The initial dataset was splitted randomly in three groups of equal size (k = 3). Two groups were used to calibrate the model and the third group was used as a test set to evaluate predictions. The concordance between predictions and observations was measured with the Spearman rank correlations (Rho) ranging from −1 to +1 [41]. The correlation is considered to be weak for −1 < Rho < 0.25, fair for 0.25 < Rho < 0.50, moderate for 0.50 < Rho < 0.75 and strong for Rho > 0.75 [42]. Cross-validation was repeated 100 times per model to obtain a variance of Rho. The predictive performances between the models were compared using Wilcoxon rank-sum non-parametric tests.
In order to assess whether the combination of multiple NDVI time periods lead to higher predictive performances than the use of a single date, a consensus approach was also followed. Because NDVI variables of different time periods were often highly correlated (Rho > 0.70) they could not be included simultaneously into a single model [33]. Alternatively, the consensus approach provides means to combine ensembles of model predictions, and to test if the predictive performance is higher after combination. This approach is usually adopted in species distribution modelling to combine outputs of single modelling techniques and to reduce uncertainty of predictions [43]. Here, the outputs of single-date models were combined with the Weighted Average (WA) consensus method [43]. The predictions were calculated as the mean of predictions of single-date models weighted by their individual predictive performance (mean Rho). The concordance between predictions and observations was evaluated in the same way than the models based on single date.
Spatial autocorrelation of residuals was tested for each best model by computing non parametric spline correlograms using the ncf package version 1.1-3 for R [44]. Neighbours were considered until a distance of 300 km. Because a positive inherent autocorrelation was observed only for the very close surveyed plot squares (Spatial autocorrelation < 0.2 within 20 km; average distance of the closest square = 8.15 ± 7.01 km; Appendix B), we considered the potential effect as not significant in the analysis [45].

Comparison of Single-Date NDVI Models with DHI According to Species Richness Metrics
Bird-habitat models showed different explanatory performances according to the species richness metrics, data types and acquisition periods considered ( Table 2). Whatever the time period of NDVI, climate always had a significant effect on the four species richness metrics. For the models based on climate only, the explanatory power is high for woodland species (D 2 = 23%) but more limited for the other groups, with the lowest explanatory power for farmland species (D 2 = 6%). Adding NDVI or DHI variables to climate increased the performance of the models. The best fits were observed for woodland (D 2 = 45% at 250 m) and farmland (D 2 = 35% at 250 m) birds. The lowest values were found for generalist species (from D 2 = 7% to D 2 = 14% at 250 m). The DHI-based models showed similar performances in comparison with the best NDVI-based models for woodland birds (D 2 = 43%), generalist species (D 2 = 13%) and total species richness (D 2 = 20%). By contrast, the explanatory power was lower for farmland species (D 2 = 21%). The same conclusions between the DHI and NDVI-based models were obtained whatever the spatial resolution of NDVI. Regarding NDVI variables, average NDVI was more important than standard-deviation NDVI for explaining species richness metrics. For DHI, average DHI cv (i.e., seasonality, based on the coefficient of variation of annual fPAR) and average DHI cum (i.e., productivity, based on the cumulative annual fPAR) contributed most (see Appendix C).
Considering the effect of NDVI seasonal variation on the model performances, whatever the landscape context, late winter (NDVI 6-Mar) was the least accurate time period for every species richness variable. The optimal NDVI time period varied according to the bird species richness metric ( Table 2). For woodland species richness, end summer (NDVI 14-Sep) was the best time period, but only slight difference in accuracy was observed with the other periods, except for winter (NDVI 6-Mar). For farmland species richness, late spring (NDVI 26-Jun) was the most adapted. For generalist and total species richness, several periods were similarly accurate. The seasonal variation of NDVI was a lower impact on these response variables since the contribution of NDVI was limited (Appendix C).
The response shapes to the NDVI variables were different according to the bird group considered (Figure 4). For woodland birds, species richness increased continuously with the average NDVI variable. By contrast, farmland species richness was quite stable for low values of average NDVI but decreased when average NDVI reached 0.7. Total richness increased until the same threshold (average NDVI of 0.7) and decreased after this value. The curves were constructed after selecting the best NDVI time period for each bird group. In terms of predictive performance, the best predictions were obtained for woodland and farmland birds (maximum mean Rho = 0.70 and 0.53, respectively) ( Figure 5). Generalist and total species richness showed fairly good but lower predictive accuracies (maximum mean Rho = 0.36 and Rho = 0.39, respectively). Comparing to the single-date models, the consensus NDVI approach showed superior predictive performances only for generalist species richness. The predictive accuracy of DHI was slightly but significantly lower than the best NDVI single-date model for all the bird groups. The national spatial distributions of woodland and farmland species richness were predicted using their respective best model, showing a consistent pattern over France. Where the farmland bird species richness is high, the richness is low for woodland birds (Figure 6). In terms of absolute value, caution must be taken because of projection uncertainty in some areas (limited bird plots in mountain and Mediterranean regions used for model calibrations) [46].

Comparison of NDVI Single-Date Models According to the Land Cover Context
The effect of seasonal changes of NDVI on models performance depends on the bird community but is also modulated by the dominant landscape composition existing in the bird plots. In other words, for the same group of birds, the best NDVI acquisition date may vary according to the dominant land cover. This effect was evaluated for woodland and farmland birds only.
In forested landscapes, with plots dominated by woody vegetation, late winter (NDVI 6-Mar) was the least accurate time period to explain woodland species richness (∆D 2 = 5% in comparison to the best date; Table 3). The other periods fitted better with close explanatory power. For farmland species richness, late spring (NDVI 26-Jun) or late summer (NDVI 29-Aug) were the most appropriated.
In agricultural landscapes, with plots dominated by winter crops or summer crops, end summer (NDVI 14-Sep) and begin of autumn (NDVI 30-Sep) explained better farmland and woodland species richness respectively (Table 3).
In grassy landscapes, with plots dominated by permanent grasslands, late summer (NDVI 29-Aug) fitted better for woodland birds. For farmland birds, begin of autumn (NDVI 30-Sep) remained the most adapted.
In a general way, landscape composition in the bird plot had an impact on the best NDVI time period for explaining woodland and farmland birds. However, end summer (NDVI 14-Sep) and begin of autumn (NDVI 30-Sep) appeared to be the most appropriate in majority (Table 3). Predictive accuracies showed close results with explanatory performances (Figure 7). Difference appeared only for farmland birds in agricultural landscapes composed of winter crops but Rho values were generally low for this case. Comparing to the single-date models, the consensus approach enabled to improve the predictions significantly for only two cases (SR woodland in grassy landscapes and SR farmland in agricultural landscapes with summer crops) (Figure 7).

Single-Date NDVI Is as Accurate as DHI or Multi-Date NDVI to Explain and Predict Bird Species Richness
Our study is the first to evaluate the potential added value of DHI components, compared to NDVI, to model the relationships between bird community and landscape structure. The results suggest that the integrative DHI components which summarizes vegetation productivity and heterogeneity over the year are not the best predictors to explain and predict diversity patterns of breeding birds over France. We found very close performances between NDVI and DHI-based models for woodland birds, generalists and total species richness. However, we observed that single-date NDVI models fitted and predicted better than DHI model for farmland birds (with a difference around 14% for the explained deviance). This is true whatever the spatial resolution of the NDVI product. NDVI at 500-m was highly correlated with NDVI at 250-m (r > 0.92 between average NDVI variables for all dates). Contrary to our assumption [21], difference in spatial resolution had no effect on models' performance.
At the scale of analysis (bird plot of 4 km 2 ), average and standard deviation of NDVI 250-m and NDVI 500-m are similar.
A possible explanation of the non-superiority of the DHI could be the difference in proxies of vegetation productivity. The DHI is based on the fPAR (derived from a physically based model) and not on NDVI which may saturate at high biomass. But a recent study showed that DHI derived from vegetation index (such as NDVI or EVI) was weaker to explain bird species richness than DHI computed from MODIS fPAR, LAI or Gross Primary Productivity (GPP) [29].
Previous studies testing the predictive power of DHI on bird species diversity were based on very large sample units such as USA ecoregions (around 90,000 km 2 in average) that are homogeneous in terms of environmental conditions [29,48]. Compared to our analysis, similar degrees of correlation were observed at this spatial extent with variations depending on avian functional groups. In [29], the multiple regression model predicted species richness of woodland birds using the fPAR-derived DHI with an adjusted coefficient of determination (%R 2 adj ) equals to 0.58. In [48] the explained variance was weaker (%R 2 adj = 0.37). In both cases, cumulative DHI was most predictive which can be referred to the species-energy theory [49]. In our study, the sample unit of investigation is substantially smaller (4 km 2 ) and match with the bird survey. At this fine scale, birds mainly respond to local landscape factors such as resource availability, habitat amount and habitat diversity [50]. The results show that DHI discriminates the habitat composition and heterogeneity less accurately at this scale than an appropriate single-date NDVI. This is especially true in agricultural landscapes where differences in DHI between crops, grasslands and woodlands may be less pronounced than NDVI at best date because of the integrated greeness over the year. Therefore, the average and standard deviation NDVI variables (proxies of dominant habitat identity and spatial heterogeneity) explain more bird community variance than cumulative and variation DHI (proxies of annual productivity and seasonnality). This finding is confirmed by the fact that the consensus approach that combine multiple time periods of NDVI has no real effect on models performance. This is also in line with previous works which demonstrated that NDVI-based measure of seasonal vegetation productivity is a stronger predictor than annual NDVI [17,51].

Best Single-Date NDVI Is Not Synchronous with the Acquisition Time Period of Bird Data
In most ecological studies, it is usual to make the acquisition date of images consistent with the ground observations in order to avoid misinterpretations [14,17,52,53]. However, as showed by the results of this study, when NDVI is used as proxy to estimate species diversity at fine scale, this common sens rule must be revised since NDVI is time-dependent.
This observation supports the conclusions of a previous study carried out in a small area of southwestern France [21]. The best NDVI time period to predict bird diversity patterns at the survey level is the one that enable to separate the mixture of crops, open grasslands and ligneous habitats. These components explain bird requirements for breeding and foraging according to functional groups. Depending on the season, differences between habitats are more or less apparent in the NDVI images ( Figure 8). Because NDVI in winter has lower values in deciduous forests, this time period leads to underestimate woodland species richness. For other seasons, models performance is equivalent. For farmland species, overestimation appears for time periods where semi-natural habitats are confused with crops because of close values of NDVI. This is the case for winter and summer crops before senescence or harvesting (at the end of june and august, varying by the location). After soil tillage, cultivated areas and other habitats are better discriminated.
Until now, this time effect has paid little attention in species distribution models based on continuous remotely sensed variables. Several studies evaluated if seasonal bird richness (breeding vs. wintering) was linked to change in productivity during the year [17,53] but few of them analyzed the optimal matching between the seasonal NDVI and the bird survey although asynchrony was already observed [22,52]. This is probably due to the macro ecological hypothesis that underpins these studies and which supposes that productivity determines species richness [54]. In our case, average NDVI reflects dominant habitat identity in each plot and this dominance is best displayed in late spring or begin of autumn. This variable was especially relevant for explaining woodland species richness and farmland species richness (Appendix D). Woodland species richness, which is mainly composed of tree-nesting species feeding on invertebrates (e.g., Certia brachydactyla, Erithacus rubecula, Sitta europea) increased constantly with NDVI. Farmland species richness was constant until a NDVI value of 0.7 and then decreased. The lowest values of NDVI correspond to the most open landscapes (Appendix E) that facilitate farmland ground-nesting species avoiding high vegetation structures (e.g., Alauda arvensis, Alectoris rufa, Coturnix coturnix). Intermediate values of NDVI represent heterogeneous landscapes with grasslands and hedgerows that favor farmland species nesting in hedgerows and requiring complementation of habitats (e.g., Carduelis cannabina, Emberiza cirlus, Emberiza citrinella) [6]. Landscape heterogeneity, through standard deviation NDVI, also explained a part of the species richness of the two bird groups but to a lesser extent. Heterogeneity may increase differentiation of ecological niches allowing different species to coexist in a same landscape [14,55,56]. However, the effects of heterogeneity are related to the degree of habitat specialization for the species. Highly specialized bird species tend to be negatively affected by landscape heterogeneity [57]. A more refined discrimination between bird species of each community, according to their habitat specialization, may have led to different conclusions [58].

Best Single-Date NDVI Is Context-Dependent, Varying According to the Dominant Land Use
Our results showed that best single-date NDVI for modeling bird distribution patterns is time-dependent but is also modulated by the landscape context of the bird plots. In landscapes dominated by winter crops or summer crops, we found that end summer (14 September ) and begin of autumn (30 September) were the best periods to explain farmland and woodland bird species richness over France. For other seasons, crop cover with high NDVI values affect negatively the estimation of the species richness (leading to an under or overestimation depending on the bird group considered). In landscapes dominated by grasslands, late spring (26 June) and late summer (29 August) were the most adequate periods for SR woodland. This is also the case in woody landscape to explain SR farmland. These time periods match with the peak of woodland and grassland vegetation productivity leading to better accurate estimations.
Because the proportion of the dominant land use in the surveyed square was sometimes low in the analysis (although a majority proportion was ensured in bird plots), a deeper analysis should be further developed to confirm the results, using a most exhaustive dataset.

Adding Habitat to Climate Variables Improve Predictions of Local Bird Patterns at National Scale
We observed that adding NDVI or DHI data reflecting local habitats in bioclimatic models greatly improve the performances of species richness estimations. The utility of habitat variables was particularly important for woodland and farmland specialists, and much less for generalist species whose individuals can find their resources in a wide range of habitats. This result is in line with previous studies showing the role of climate as driver of bird species distribution at macro-scale, combined with land cover at finer spatial resolutions [59,60]. However, the role of climate as causal determinant of species' geographic distributions is still in debate [61,62]. Positive correlation with climate could appeared for non-climatic reasons, due to the existence of spatially structured processes which are themselves not necessary related to coarse environmental variables [63]. At fine resolution, landscape and habitat composition which is highly variable between plots effectively increases the performance of species richness models. NDVI data can thus be of great value for building bird models, especially for specialist species whose populations decline drastically [64].

Conclusions
As confirmed in this study, continuous (unclassified) predictors derived from satellite imagery contribute to build efficient SDM at large spatial scale [7]. However, because of temporal variability of remotely sensed predictors over the year, acquisition time periods of the data affect the models' performance. Here, we found that: • Compared to a combination of multiple NDVI time periods or to the components of DHI, performance of bird models is higher with single time period of NDVI variables. • Models' performance is time and context-dependent. At the French scale, end summer (NDVI 14-Sep) and begin of autumn (NDVI 30-Sep) are in most cases, the best time periods to explain and predict woodland and farmland birds species richness. • Best time period of NDVI is the one that better reflect dominant habitats in the bird survey plots.
• Habitat components at finer scale have to be combined with climate variables at macro scale to explain and predict patterns of habitat specialist species.
These findings suggest that the most powerful approach is the simplest one. Only one image of appropriate time period is required for mapping bird species richness. This makes the use of unclassified remotely sensed data very operational for bird monitoring programs. Further work is needed to confirm the conclusions from multiple years and to verify the stability of the relationships. Predictions could also be improved using the new Sentinel-2 time series of higher spatial resolution.