Comparing Global Sentinel-2 Land Cover Maps for Regional Species Distribution Modeling

: Mapping the spatial and temporal dynamics of species distributions is necessary for biodiversity conservation land-use planning decisions. Recent advances in remote sensing and machine learning have allowed for high-resolution species distribution modeling that can inform landscape-level decision-making. Here we compare the performance of three popular Sentinel-2 (10-m) land cover maps, including dynamic world (DW), European land cover (ELC10)


Introduction
The Anthropocene has heralded an unprecedented loss of biodiversity, primarily due to changes in land use and land cover, climate change, pollution, (over-)exploitation, and biological invasions [1]. In response, governments have established frameworks that address biodiversity loss, including the United Nations' (UN) Sustainable Development Goals (SDGs, 2030 Agenda) as well as the Aichi biodiversity targets (Strategic Plan for Biodiversity 2011-2020) and the Post-2020 Global Biodiversity Framework of the Convention on Biological Diversity. Recently, the UN established the statistical standards for ecosystem accounting (EA) under the System of Environmental-Economic Accounting (SEEA), which require countries to account for changes in ecosystems over time [2]. To achieve the SDGs, meet biodiversity conservation targets, and to account for ecosystem changes, we require monitoring and evaluation tools that are both globally available and locally relevant [3].

Methods
To compare the utility of Sentinel-based land cover maps in species distribution modeling, we chose wild bee species as a model system. We do this for three main reasons, which are elaborated on below: (1) bees are keystone species in grassy ecosystems globally that are good indicators for ecosystem condition; (2) they are experiencing significant local declines in population numbers; and (3) they have a limited home range and are dependent on local resources for survival, so we expect their distribution to respond strongly to landscape-scale land use gradients in the Sentinel-based land cover maps.
Insects in general and wild bees in particular are key components of many ecosystems and provide important contributions to people [19,20]. For example, the economic value of pollinating insects is considerable, estimated at EUR 153 billion globally [21]. Although pollination is often associated with domestic insects such as the honey bee (Apis mellifera), a large and diverse community of wild bee species contributes significantly to the pollination of wild plants and crops [22,23]. Insect abundance, biomass [24], and diversity [25] are declining in some regions due to urbanization, deforestation, climate change, pesticide use, Remote Sens. 2023, 15,1749 3 of 13 and invasive species [26]. The same is true for bee species [27]. Although the diversity of wild bees is predominantly driven by temperature at global and regional scales, habitat and landscape characteristics are important drivers at local scales [28,29]. Bees are central place foragers that travel back and forth from a nesting site to collect resources [30], and therefore, the abundance and diversity of wild bees correlate to the diversity and abundance of floral resources (vascular plants) in the landscape [28,31]. Land cover types such as semi-natural grasslands, including hay and flower meadows, provide ample floral [32] and non-floral resources to wild bees [33]. Non-floral resources include nesting sites, and although some wild bee species nest in dead wood, the majority use sandy soil sediments as nest sites [34]. Consequently, wild bee diversity can be predicted by the availability of suitable land cover types for resources (i.e., semi-natural grasslands) and for nests (i.e., sandy soils), in combination with climatic variables [29].

Solitary Bee Surveys
We surveyed solitary bee communities in 72 traditionally managed (mowed) seminatural grasslands distributed along a climatic gradient from south-eastern Norway to mid-Norway ( Figure 1A). In 2019, we sampled 32 semi-natural grasslands in south-eastern Norway [11], adding another 20 study sites in 2020 [29]. To capture potential influences of climatic conditions on solitary bee diversity at regional scales, solitary bee communities were sampled in another set of 20 semi-natural grasslands in mid-Norway in 2021. All surveys were conducted using pan traps, which are an efficient method for surveying wild bee communities [35], in particular when the aim is to survey the solitary as opposed to social bee species [36]. Each pan trap consisted of three white plastic soup bowls, coated with fluorescent yellow or blue, or left white, mounted on a fence pole at the height of the surrounding vegetation [11,29]. We deployed 3 traps per site in 2019 and 2 traps per site for sites sampled in 2020 and 2021. The number of traps per site was reduced from 3 to 2 after 2019 to allow more sites to be sampled. This resulted in 176 samples (the sum of traps across all sites) distributed across 72 study sites. Within sites, traps were always placed at least 20 m apart to avoid inter-trap competition [37].

Land Cover Maps
Data pre-processing and extraction of land cover maps took place in the Google Earth Engine [38]. The AR5 map obtained from the Norwegian Institute of Bioeconomy Research [18] is provided as a vector map at 1:5000 scale, which we rasterized at 5 m resolu- In all years, sites were sampled in May, June, and July at similar times of day and on days with weather conditions that are optimal for bee activity (i.e., low wind and temperatures above 15 • C). For each sampling period, we activated the pan traps at a site by mounting the bowls and filling them with water and a drop of detergent. Sampled bee specimens were collected after 48 h. In 2019, sampling was initiated on 13 May, 21 June, 9 July, and 23 July [11]. In 2020, sampling was initiated on 13 May, 25 May, 14 June, and 16 July [29]. In 2021, sampling began on 29 May, 22 June, and 24 July. Collected bees were stored in 96% laboratory ethanol prior to pinning and identification. Voucher specimens are stored in the entomological collections at the Norwegian Institute for Nature Research. We tallied the number of solitary bee species sampled per trap across the season.

Land Cover Maps
Data pre-processing and extraction of land cover maps took place in the Google Earth Engine [38]. The AR5 map obtained from the Norwegian Institute of Bioeconomy Research [18] is provided as a vector map at 1:5000 scale, which we rasterized at 5 m resolution. AR5 is updated every 5-8 years and therefore represents a mosaic of years across the country. WC and ELC10 are available in Google Earth Engine for the years 2020 and 2018, respectively. However, DW is provided as a collection of classified Sentinel-2 images with less than 35% cloud cover. To generate an annual land cover composite for 2020 comparable with WC and ELC10, we calculated the mode predicted land cover class in the image band named "label" across all DW images during June, July, and August. For all land cover maps, we used a focal mean function to calculate the proportion of grassland habitat within 250 m of each pixel in the study area ( Figure 1) at 10 m resolution. We used the 250-m radius because solitary bee diversity has been shown to respond strongly to habitat availability at this spatial scale [39]. To isolate grassland pixels, we used the "grassland" class from WC and ELC10 and the "grass" class from DW, which are both defined as areas dominated by natural herbaceous vegetation, including grasslands, prairies, steppes, savannahs, and pastures. For AR5, we used the "innmarksbeite" and "åpent fastmark" classes, which are defined as open ecosystems dominated by herbaceous vegetation and often used for extensive grazing [18]. We used a radius of 250 m as solitary bee species richness has previously been shown to respond to habitat area at this spatial scale [39].

Modeling
Data modeling and visualization were performed in R [40]. As in [29], we included predictor variables related to climatic conditions, habitat availability, and distances to high-quality nesting substrates for below ground nesting bees, which account for the vast majority of solitary bees. In [29], the spatial variation in climatic conditions was estimated using a digital elevation model (DEM) together with latitude. Here, we used the average temperature for the warmest quarter (June, July, and August), during the current 30-year climate reference period (1990-2021), calculated from daily estimates and interpolated station measurements at a 1 km resolution from the Norwegian Meteorological Institute's database [41]. We used the average temperature during the warmest quarter instead of the annual mean temperature because high winter temperatures along the coast result in annual mean temperatures not reflecting the gradient in ambient temperature experienced by bees during their main activity periods in Norway (spring to autumn). Using modeled climate data instead of DEMs to estimate the effects of temperature on bee diversity enables one to project changes in bee diversity as a function of future climate scenarios. In [29] we used a habitat suitability model to estimate the potential habitat area surrounding each pan trap at a 60-m radius. In contrast, in the current study, we use the proportion of pixels identified as "grassland" by satellite-derived land cover maps as estimates of habitat availability within a 250-m radius surrounding each trap. A benefit of using satellitederived grassland classifications instead of estimates of habitat suitability maps that may partly rely on existing land cover products is that one reduces the number of modeling steps required to produce maps of pollinator diversity from remote sensing data. In addition, we Remote Sens. 2023, 15, 1749 5 of 13 used the geographic distance to areas with soils mapped as having a high water infiltration capacity by the Norwegian geological survey [42], as such areas are typically located on sandy soils, which is the preferred nesting substrate for the majority of Norwegian soil nesting bees.
We ran and compared eight models of solitary bee species richness that all followed the general formula: Solitary bee species richness (survey year + average temperature during the warmest quarter + Distance to sandy soils + proportion of grassland within a 250-m radius).
The eight models differed in terms of the data source (map) used to estimate the proportion of grassland around each pan trap and the type of model used to predict solitary bee species richness. The survey year was included as a categorical variable to account for inter-annual variation in bee species richness as well as for annual differences in climatic conditions, which could influence the number of species sampled in a trap. For each land cover map, we trained a Random Forest (RF) regression model [43] and a Bayesian Regularized Neural Network (BRNN) model [44] in Caret [45] in R. BRNN models were tuned by selecting the number of neurons (1, 2, or 3) that resulted in the lowest root mean square error (RMSE) following 25 bootstrap resamples of the training data. RF models were tuned by selecting the mtry (number of parameters tested at each node) and split-rule variance or extra trees that resulted in the lowest RMSE, following 25 bootstrap resamples. We adopted leave-one-out cross-validation (LOOCV) to assess model predictive performance due to the small number of study sites. Small training datasets can result in large variances in model performance when using traditional training-testing splits (e.g., 70% training and 30% testing) before model fitting, compared to LOOCV [46]. In LOOCV, each model was iteratively trained and tuned on data from 71 study sites and then used to predict solitary bee species richness in pan traps from the one held-out site. Predictive performance is calculated across all LOOCV iterations. Following this nested cross-validation procedure ensures independence between the data used for tuning the models and the data used for final validations.
For each land cover-map model prediction, we evaluated the predictive power in terms of the Pearson correlation coefficient (Cor), the root-mean-square deviation (RMSE), and the mean absolute error (MAE) between observed values of solitary bee species richness and the solitary bee species richness predicted by the model. As an alternative form of model validation, we also tested the ability of model predictions of bee species richness to predict the variance in occurrence of solidary bee records obtained from citizen science data. We first downloaded all post-2015 solitary bee species observations from the Global Information Biodiversity Facility (GBIF) that intersected our study area (xmin: 9.836, xmax: 12.719, ymin: 58.909, ymax: 63.962) and had a GPS error of less than 50 m (n = 2111). To reduce spatial bias in the data-i.e., that some areas have been surveyed more frequently or intensively than others-we only included one occurrence per square kilometer. We randomly sampled 10,000 pseudo-absences from within our study area and used binomial generalized linear models to quantify how well the predicted species richness scores explained the variance in the GBIF presence-absence data. We used the Akaike information criterion (AIC) to compare the explanatory capacity of the different models (δAIC > 2). Finally, we used partial dependency plots through the R-package pdp [47] as a means to visualize the estimated marginal effect of the proportion of grassland, using the different land cover maps, on solitary bee species richness.
In summary, our modeling workflow included 16 unique models for combinations of reference data (survey, GBIF), model type (RF, BRNN), and land cover dataset (AR5, DW, ELC10, WC). Each model was iterated 100 times in order to estimate the mean and standard deviation in model performance metrics.

Results
We did not find that models using the proportion of grassland estimated from the vector-based Norwegian land cover map (AR5) outperformed models where grassland had been estimated from satellite-derived land use models with 10 m resolution. Specifically, satellite-derived maps exhibited an average RMSE of 2.87 ± 0.03 (±standard deviation), whereas the Norwegian AR5 map produced models with a RMSE of 3.02 ± 0.02 (Figure 2). Similarly, for the AR5-based BRNN models, the correlation coefficient between observed and predicted solitary bee species richness (SR) was slightly lower (Figure 2A) than for the satellite-based models ( Figure 2B-D). This was also the case for the Random Forest (RF) models ( Figure 2E-H).  Among the satellite-derived models ( Figure 2B-D,F-H), there were only marginal differences in their performances as predictors of solitary bee SR. Models using grassland habitat from DW performed best (RMSE = 2.8 ± 0.03; averaged across RF and BRNN models), followed by ELC10 (RMSE = 2.85 ± 0.03) and WC (RMSE = 2.87 ± 0.02). Grassland habitat was ranked the third most important predictor variable in BRNN and RF models ( Figure 3). The mean summer temperature and sampling year were the most important predictors in the models. All models, independent of data source or model type, produced partial dependence plots that corroborated a positive association between grassland habitat and solitary bee SR (Figure 4). The association was less distinct in the AR5 model ( Figure 4A   When validating the model predictions of bee SR against citizen science data on solitary bee occurrences using generalized linear models ( Figure 5), we found that the order of performance was changed. We found that ELC10 performed best (AIC = 2278 ± 4), followed by WC (AIC = 2367 ± 3), and DW (AIC = 2376 ± 3). The BRNN models performed slightly better than RF models, both in terms of leave-one-out cross-validation ( Figure 2) and in terms of explaining the occurrence of solitary bees from GBIF ( Figure 5). Therefore, we used the BRNN models to generate wall-to-wall prediction maps of bee SR over the study region ( Figure 6). A qualitative visual comparison shows that the broad-scale spatial patterns of bee SR predictions are similar between models, with bee SR increasing with north-south temperature gradients and along populated valleys where extensive grazing practices have established semi-natural grassland patches. At the landscape scale (Figure 7), AR5 predictions show less spatial variation than the satellite-derived maps. All models appear to pick up the grassland habitat adjacent to the runways at Gardermoen International Airport; however, ELC10 appears to pick up the most habitat in the agricultural landscapes southwest of the airport ( Figure 7C).     neural network (BRNN) and random forest (RF) models. Effects plots are shown along with AIC scores from binomial generalized linear models of solitary bee species occurrence data from GBIF as a function of predicted bee SR from models trained with land cover data from (A,E) a Norwegian map (AR5), (B,F) dynamic world (DW), (C,G) European land cover (ELC10), and (D,H) world cover (WC). Points in the figures show the mean occurrence of solitary bee records calculated within bins of one decimal. Individual GLMs were run for each of the 100 spatially filtered datasets of solitary bee records in order to calculate the average AIC and its standard deviation.

Discussion
Species distribution modeling that incorporates high-resolution satellite data is still in its infancy [13], yet the application to solitary bee species richness presented in this study confirms its potential. Our results show that the Sentinel satellite-based land cover maps outperformed a regional manually digitized land cover map over southern Norway (AR5) in predicting solitary bee species richness (Figure 2). While the use of satellite imagery to map vegetation types or individual forest species is common [48], using derived

Discussion
Species distribution modeling that incorporates high-resolution satellite data is still in its infancy [13], yet the application to solitary bee species richness presented in this study confirms its potential. Our results show that the Sentinel satellite-based land cover maps outperformed a regional manually digitized land cover map over southern Norway (AR5) in predicting solitary bee species richness (Figure 2). While the use of satellite imagery to map vegetation types or individual forest species is common [48], using derived products to predict habitat suitability for animal or plant species that are not directly visible in satellite images is far less common [11,49]. Here we show that the availability of grassland habitat within 250 m, measured from satellite data, is positively associated with solitary bee species richness ( Figure 4) and is therefore predictive of solitary bee occurrences at regional scales ( Figure 5).
We found larger differences in the predictive capacity of grassland habitat derived from manually mapped versus satellite-derived land cover maps than differences between the satellite-based maps themselves (Figure 2). The AR5 map uses a minimum mapping unit of 2000 m 2 , which results in very small grassland fragments being subsumed into a broader land cover class [18]. For example, road verges or small urban parks will be classified as "built" in the AR5 maps. However, such small grassland patches can harbor significant floral resources for bees, and the fact that AR5 does not map these areas may therefore be why AR5 was less predictive of species richness than satellite-based maps. Furthermore, AR5 is not as up-to-date as the satellite-based maps and may misrepresent the conditions on the ground during 2019, 2020, and 2021, when the field surveys were conducted. In contrast, ELC10 and WC use a minimum mapping unit of 100 m 2 , while DW uses 2500 m 2 . At the landscape scale, it is evident that predictions of bee species richness with ELC10 were more spatially heterogeneous than AR5 or DW (Figure 6), probably due to its smaller minimum mapping unit and ability to detect smaller grassland patches. This may have contributed to its greater predictive capacity compared to AR5. This also explains why, when validating our models using citizen science data from GBIF, ELC10 outperformed DW. In contrast to the survey dataset, GBIF data are spatially clustered and biased toward urban landscapes, which are easily accessible but also have complex landscapes. Due to ELC10's small minimum mapping unit, it captures the landscape complexity more than DW does and is able to predict the GBIF bee SR better.
The accuracies of the solitary bee species richness models presented here are arguably high enough to make data-driven management decisions at the landscape scale. The average RMSE of 2.87 means that one can at least distinguish very species-rich areas (maximum species richness of 16 in our study) from species-poor areas (zero species). The differences in model accuracy between satellite-based grassland maps were marginal (RMSE difference of 0.04); however, when visualized at a landscape scale (Figure 6), small nuances may have important implications for management and decision-making. For instance, a ranking or prioritization of grassland patches for receiving conservation subsidies based on the maps in Figure 6 may yield different results depending on the data source used. Therefore, post-stratified accuracy assessment of species distribution models in specific landscapes may be necessary before they can be adopted in practice [50].
Based on several limitations identified in our study, we outline avenues for further research on the integration of high-resolution satellite data in species distribution modeling. Firstly, it is not necessary to use derived products such as land cover maps, as we have carried out here. Instead, one can use the spectral signatures themselves from a satellite image to calibrate distribution models [13], although this would not allow ecologists or policymakers to relate species-rich areas to specific land use types. Secondly, maps with a higher thematic resolution than those used in this study would produce more detailed species distribution maps [51]. For example, all four maps tested here contained a single broad category for grassland, without distinguishing between intensively managed grasslands and extensively managed grasslands, such as the mowed meadows from which we sampled solitary bees [11,29]. Therefore, measuring aspects of ecosystem condition or use, such as grassland use intensity [52], might further improve the accuracy of species distribution models. Thirdly, we did not explore how accurately satellite-based prediction models can detect real changes in species distributions over time because we did not implement a field survey design that was comprehensive enough to capture changes in species ranges over time. However, we know from earlier studies that changes in land cover that are detectable from space are direct drivers of species range shifts [53]. To this end, DW is probably the most suitable Sentinel-based land cover dataset to quantify land cover and use dynamics [54] because of its continuous updates and delivery and would be well-suited to such dynamic species distribution modeling. This also strengthens the call for investment in long-term biodiversity monitoring programs so that satellite-based distribution models can be calibrated and validated with in situ data [55].

Conclusions
The proliferation of high-resolution earth observation data and derived land cover products provides scope for mapping biodiversity distributions with models that are both locally relevant for decision making and scalable to the globe. Here we found that globally available Sentinel-based land cover maps can improve upon manually digitized regional land cover maps for predicting the richness of solitary bee species in southern Norway. The differences in predictive performance between DW, WC, and ELC10 were marginal; however, at the landscape scale, the smaller minimum mapping units of WC and ELC10 allow them to resolve smaller habitat patches, which are reflected in the landscape variations in predicted species richness. Furthermore, the rich time series provided by maps such as DW (from 2015 to present) offer unique opportunities to model short-term changes in species distributions in response to land use changes if paired with in-situ temporal monitoring data. We conclude that the use of satellite-derived land cover maps can facilitate high-resolution species distribution models that can guide decision-making relevant to landscape ecology. To this end, future modeling efforts should be aimed at those species that perform key roles in ecosystems, are indicators of ecosystem status, and support nature's contribution to people.