The Role of Remote Sensing Data in Habitat Suitability and Connectivity Modeling: Insights from the Cantabrian Brown Bear

Abstract

Ecological modeling requires sufficient spatial resolution and a careful selection of environmental variables to achieve good predictive performance. Although national and international administrations offer fine-scale environmental data, they usually have limited spatial coverage (country or continent). Alternatively, optical and radar satellite imagery is available with high resolutions, global coverage and frequent revisit intervals. Here, we compared the performance of ecological models trained with free satellite data with models fitted using regionally restricted spatial datasets. We developed brown bear habitat suitability and connectivity models from three datasets with different spatial coverage and accessibility. These datasets comprised (1) a Sentinel-1 and 2 land cover map (global coverage); (2) pan-European vegetation and land cover layers (continental coverage); and (3) LiDAR data and the Forest Map of Spain (national coverage). Results show that Sentinel imagery and pan-European datasets are powerful sources to estimate vegetation variables for habitat and connectivity modeling. However, Sentinel data could be limited for understanding precise habitat–species associations if the derived discrete variables do not distinguish a wide range of vegetation types. Therefore, more effort should be taken to improving the thematic resolution of satellite-derived vegetation variables. Our findings support the application of ecological modeling worldwide and can help select spatial datasets according to their coverage and resolution for habitat suitability and connectivity modeling.

1. Introduction

Biodiversity loss is among the biggest challenges that nature conservation is currently facing. As a result, many conservation efforts are being continuously made by scientists, government administrations and conservationists. Computer-based approaches are acknowledged as one of the most useful and unbiased instruments to guide the design of nature management plans and policies, with habitat suitability and connectivity models being particularly helpful to produce spatially explicit information for conservation [1,2,3,4,5,6]. Nevertheless, sufficient spatial resolution of data and a careful selection of environmental variables are essential to ensure reliable model performance [7]. Accordingly, considerable effort has been taken to increase spatial resolution. The use of remote sensing imagery is gaining relevance because of its potential high spatial resolution and its broad temporal and spatial coverage and availability [8,9,10].

Some studies have investigated the effects of using different spatial data sources to model species–habitat relations [11,12,13]. They found similar results regardless of the data sources used in the models, although slight variations in their performance may depend on the spatial and thematic (diversity of categorical classes in a spatial layer) resolution of the datasets. Other studies also support this idea by indicating enhanced performance when high-resolution remotely sensed data are used for ecological modeling [14,15]. For instance, the utilization of light detection and ranging (LiDAR) can improve models’ performance by informing about the structure of the vegetation and habitats [16,17,18]. Similarly, the use of high-spatial resolution optical imagery also yields reliable model fits [15]. However, these spatial datasets may be expensive to obtain—e.g., Quickbird, Pleiades, SPOT, Worldview, IKONOS—or only available at a national scale—e.g., LiDAR is only freely available in some countries such as Spain, Denmark or Great Britain. As an alternative, some Earth observation missions, such as Sentinel and Landsat, currently offer global open source high-resolution satellite data, including optical and radar imagery. In particular, the Sentinel-1 and 2 constellations have the advantage of having global coverage and a revisit time of six and five days at the equator, respectively [19,20,21]. These properties allow nature managers to access spatial data of their region of interest and estimate satellite remote sensing-based biodiversity variables [22] and characterize land cover regardless of the location. Additionally, managers can regularly update their results, which is beneficial for monitoring conservation measures [23,24]. Other alternatives include the pan-European spatial products of Copernicus Land Monitoring Services (CLMS). Although not being processed globally, this platform also offers free comprehensive high-resolution layers and land cover maps that can be a cost-effective option for projects focusing on Europe. Layers in CLMS are mostly derived from optical remotely sensed data, especially from Sentinel-2 and Landsat-8 as ancillary data. These layers include the Forest Type Product (FTY) which characterizes forest cover in Europe [25] and Corine Land Cover (CLC) which distinguishes 44 cover types through the interpretation of satellite imagery [26] (details of these and other layers can be seen in Appendix A). Despite all these alternatives, to our knowledge, little work has been conducted to compare how various spatial data sources that differ in availability and coverage affect habitat suitability and connectivity models.

Therefore, this study aimed to determine whether global and free satellite imagery can match up to costly or regionally restricted spatial data sources in ecological modeling. To achieve this, we developed habitat suitability and connectivity models for the brown bear (Ursus arctos arctos) in the Cantabrian Range (Spain) using three spatial datasets that differ in their level of spatial coverage and accessibility, i.e., global (a Sentinel-1 and 2 land cover map produced in a previous study [27]), continental (pan-European vegetation and land cover layers such as FTY [25] and CLC 2018 [26]) and national (Spanish LiDAR and Forest Map [28,29]). We compared and validated the predictive performance of the habitat models, and the dispersal routes and costs of displacement were estimated by the connectivity models. The results of this study can help to understand the implications of the choice of spatial data sources for ecological modeling, considering potential financial or spatial coverage limitations. This is especially important in the Cantabrian Range, where considerable conservation efforts have been taken to recover and reconnect the two endangered brown bear subpopulations, located in the east and west of the region [30,31].

2. Materials and Methods

A combination of species occurrence data, environmental layers and modeling techniques was used to address the objectives of the study (see Figure 1 for a workflow of the methodology).

2.1. Study Area

The study area is located in the Cantabrian Range (NW Spain) and covers 35,700 km², which includes the entire brown bear population of the Cantabrian Range. The mountainous zones have kept a considerable section of the study area with its original natural cover despite the human influence found in the region. The anthropogenic activity created a heterogeneous landscape with a mosaic of land covers where artificial and agricultural surfaces coexist with a diversity of grasslands, shrublands and forests. Natural areas are mainly composed of oak (Quercus sp.), beech (Fagus sylvatica) and chestnut (Castanea sativa) forests, heathlands (Erica sp.), shrublands (Cytisus sp.) and grasslands.

2.2. Species Occurrence Data

We used a comprehensive dataset of brown bear locations of the eastern and western subpopulations gathered from 2000 to 2010 by professional observers and rangers through direct observations and indirect evidence including excrements, footprints, beehive attacks and fur (details in [32]). This set of locations represents a time when both subpopulations were largely isolated and connectivity was identified as a major conservation concern for the species. This situation provides a strong opportunity to assess potential differences with far-reaching consequences in the definition of corridors and habitat areas where successful implementation of measures is fundamental. We unified the presence records and resampled them to 1 hectare pixel size to enhance the processing time. In total, we obtained 6207 brown bear locations to fit the habitat suitability models.

2.3. Environmental Variables and Spatial Datasets

We considered a set of six environmental variables to represent foraging resources, shelter (expressed by forest area and terrain ruggedness) and human pressure (expressed by building, highway and road density) to model the habitat suitability of the brown bear (Figure 1,

Table 1. Environmental variables for the three datasets according to their spatial coverage (global, continental or national). More details can be seen in the Materials and Methods section and Figure 1. Acronyms are Forest Type Product (FTY), Corine Land Cover 2018 (CLC18) and Forest Map of Spain (FMS).

Variable	Dataset	Data source
Foraging Resources	Global	Land Cover Map based on Sentinel-1 and 2 classification (10 m) [27]
	Continental	CLC18 land cover data [26] (vector) combined with FTY data for canopy cover (20 m) [25]
	National	FMS (vector) main tree species data and niche models for other relevant plant species, combined with canopy cover (25 m) derived from LiDAR data
Forest Area	Global	Forest class of Sentinel-1 and 2 Land Cover Map (20 m)
	Continental	FTY (20 m)
	National	LiDAR: ≥10% of canopy cover (25 m) [28]
Ruggedness	Global Continental National	Digital Elevation Model of 25 m resolution [28]
Buildings	Global Continental National	Topographic Map (vector) at 1:25,000 scale (MTN25) [28]
Highways Roads	Global Continental National	Road density (vector) derived with Open Street Maps data [33]

). Previous studies have shown this set of environmental variables to effectively predict brown bear presences in the Cantabrian Range [16,32]. Foraging resources and forest area variables were obtained from (1) a Sentinel-1 and 2-based land cover map of 2016–2017 (SLCM) [27], (2) the Forest Type Product 2015 [25] of CLMS, (3) the CORINE Land Cover 2018 map [26], (4) the Forest Map of Spain [29] (FMS; 1998–2007) and (5) LiDAR data (2009–2012) from the Spanish National Plan for Aerial Orthophotography [28]. See Appendix A for more details of these datasets. Terrain ruggedness and human pressure variables were constant in all the models and derived from the same data sources.

We combined the spatial data from the previous sources into three datasets according to their spatial coverage—i.e., global, continental or national—to produce the environmental variables used as predictors of habitat suitability. However, the extent of the analysis was constant (i.e., the Cantabrian Range) and the terms “global”, “continental” and “national” were only used to indicate dataset characteristics regarding spatial coverage, data availability and thematic, spatial and temporal resolution (see Appendix A for dataset details). Therefore, for the global level, we used the SLCM alone. Meanwhile, we combined CLC18 and FTY for the continental level, and LiDAR with FMS for the national level.

We resampled the environmental variables to 1 hectare pixel size to match species occurrence data resolution. We considered nine different spatial scales for each variable, each scale representing a particular brown bear’s point of view and reaction to the environment. This environmental gradient is of particular relevance for its ecological requirements, as the brown bear may behave differently when selecting resources from patch to landscape scales [34]. We calculated mean values using circular moving windows with nine different radii (0.25, 0.5, 1, 2, 4, 8, 16, 32 and 64 km). We standardized all variables by subtracting the mean and dividing by the standard deviation.

2.3.1. Foraging Resources

We calculated foraging resources as the abundance of plant species (trees, shrubs and herbs) weighted by the importance of each species in the diet of the brown bear (for details, see [32]). Importance was obtained from a previous report on the brown bear diet based on scat surveys [35]. Abundance was calculated in different ways depending on available information of each spatial dataset. For the national dataset, we used the information available in the FMS map, which includes tree species identities and total tree canopy cover and percentage of cover for each tree species. However, we replaced the canopy cover from FMS by the one derived from LiDAR data, which is more detailed due to its finest spatial resolution (25 m; Appendix A). We estimated the non-tree species abundance combining information from available floristic inventories, ecological niche models fitted through penalized logistic regression and based on climatic and lithological predictors (Appendix C) [36], and expert knowledge on the potential presence of specific plant species within plant communities identified by FMS. We used plant species that are basic for the brown bear diet (Appendix C). For the continental and global datasets, we used the abundances estimated in [35], adapting them to fit CLC18 and SCLM classes by averaging the original abundance values for each class. Additionally, for the continental dataset, we replaced the forest covers of CLC18 with the ones of the FTY due to its finer grain size.

2.3.2. Shelter

We used forest area and terrain ruggedness as proxies for brown bear shelter following the methods described in [16,32]. The definition of forest area varied slightly among datasets. For the global dataset, we obtained it from the forest class of the SLCM. We used the FTY layer for the continental dataset. For the national dataset, we applied a threshold of 10% of canopy cover of LiDAR data to determine forest pixels, which is the same threshold applied to the FTY (Appendix A; [25]). We calculated terrain ruggedness from a 25-m Digital Elevation Model [28], according to [37].

2.3.3. Human Pressure

To study human pressure, we used three variables that represent buildings and transport infrastructures: the density of buildings, highways and conventional roads. We obtained buildings data from a topographic map at 1:25.000 scale (MTN25 2015) [28] and transport infrastructures from Open Street Maps 2015 (downloaded from [33]). Transport infrastructures were divided into conventional roads and highways due to their different traffic volumes and physical limitations for the brown bear distribution. Both vector datasets were rasterized with 100 m pixel size and were then processed in the same way as other variables for the multi-scale analysis (see Section 2.3).

2.4. Habitat Suitability Models

First, we developed univariate models for every combination of environmental variable, scale (radii of the circular moving window; see Section 2.3) and dataset (global, continental and national) to find the operational scale, which is the scale that correlates best with the occurrence data (Figure A1, Figure A2 and Figure A3) [32,38,39]. For this purpose, we applied lrm and pentrace functions of the rms R package [40] to fit penalized logistic regressions with data of each scale as a unique predictor and assuming a monotonic response. The outcome binary variable was built merging bear occurrences with 20,000 pseudo-absences generated randomly as background cells inside the study area. We selected the model with the best value of Akaike’s information criterion (AIC; [41]). The scale of the best performing model was identified as the operational scale (i.e., the scale that best represents the habitat selection pattern of the species for each variable).

Second, we used every environmental variable at its operational scale to fit multivariate and multi-scale habitat suitability models for each of the three datasets. We calculated penalized logistic regression models using linear terms and no interactions among predictors, using lrm and prentrace functions of the rms package. To study the predictive ability of each of the three multivariate models, we relied on the area under the receiver operating characteristic curve (AUC; [42]), tested with ten-fold cross-validation. Next, we tested the difference between the AUC of the models by applying DeLong’s test [43]. The predicted suitability values of each multi-scale model were mapped to obtain three habitat suitability maps (one for each dataset).

2.5. Connectivity Analysis

To assess brown bear connectivity, we used an individual-based approach and the graph theory, which is one of the most common methods to study ecological connectivity [44,45,46]. Graphs are mathematical representations of the landscape that allow effectively studying complex connectivity scenarios [47,48]. Graphs are composed of nodes and links between nodes, which ultimately represent habitat patches and movement corridors, respectively. The links are preferably characterized by effective distances which consider the species cost of movement through the landscape matrix [49]. Thus, higher distances indicate a high probability of isolation between source and destination habitat areas. To compute effective distances, it is therefore necessary to have a resistance surface (i.e., a raster layer) whose pixel values express the cost of movement, usually derived from the land covers present in each cell [50,51]. We derived resistance surfaces from each of the three bear habitat suitability maps by assuming an inverse relation of movement resistance with habitat suitability using the following expression:

$$\mathrm{R}=1+\frac{\left(\max \left(\mathrm{p}\right)-\mathrm{p}\right)\times 99}{\max \left(\mathrm{p}\right)}$$

(1)

where R is the resistance value calculated for a given pixel, p is the value of habitat suitability in that pixel and max(p) is the maximum value of habitat suitability across all pixels in the study area.

The transformation of habitat suitability values resulted in three resistance surfaces with values ranging from 1 to 100. The lowest values represent pixels of high habitat suitability and low resistance and vice versa. For each of the three resistance surfaces, we computed distances between nodes in two ways according to two of the three connectivity methodologies of this study. Nodes consisted of an array of 316 pixels resulting from integrating the brown bear presence records into a grid with 5 km cells (a size previously used for the species [38,52]) to reduce spatial autocorrelation and avoid overlaps.

We used three methodologies to assess connectivity, the first two to characterize links between nodes and a third one to determine the links’ importance for connectivity. First, we used least-cost modeling to characterize links as the paths that imply the minimum effective distance between two nodes (i.e., least-cost paths), which are more likely to be used by animals due to lower dispersal efforts and mortality risk [53]. We used Linkage Mapper software to compute least-cost paths and their effective distances [54].

Second, we used circuit theory [55] to characterize links. Circuit theory allows the calculation of currents. Currents represent net probabilities of movement which are estimated from random walks between nodes and resistance surfaces. Corridors are characterized as a set of cells with high probabilities of being used—i.e., high current flows. As in least-cost modeling, circuit theory measures effective distances but incorporates multiple pathways connecting nodes (i.e., resistance distances). Therefore, several available dispersal routes among locations can be identified including optimal and suboptimal ones. We applied circuit theory with Circuitscape v4.0.5 [56].

For the third methodology, we used the probability of connectivity index (PC). This index is based on the habitat availability concept, which assumes that connectivity not only happens between habitat patches, but also inside them [48]. The probability of successful connection between two nodes is calculated as a decreasing exponential function of the effective distances between them. In our case, the effective distances between nodes were calculated through least-cost modeling. Then, these values are summed up to consider all the connections present in the landscape (for details, see [49]). Therefore, it is possible to calculate the contribution of each link to overall connectivity by calculating the difference in PC (dPC) when each link is removed, those links with a higher dPC being more important [57]. To calculate dPC values, we used the Link Removal function of Conefor software [58]. To estimate the probabilities, we used the median dispersal distance of the species, which represents a 0.5 probability of connection between two nodes separated at that distance. We considered a dispersal distance of 119 km for the brown bear [59], which we multiplied by the mean value of each resistance surface to obtain the effective dispersal distances.

We compared the results derived from the three resistance surfaces (global, continental and national) when applying the three methodologies (least-cost modeling, circuit theory and PC index). First, we obtained connectivity maps developed from each methodology and visually compared the distribution and location of movement pathways. Second, we compared the correlation of the connectivity metrics obtained from the three methodologies (effective and resistance distances and dPC values) between resistance surfaces. Last, we applied the three methodologies to a genetic resistance surface (genetic model) and compared the results to those obtained from the three previous resistance surfaces. The genetic resistance surface was obtained from brown bear genetic samples (17 polymorphic microsatellite loci) integrated with landscape variables following a genetic-multiplicative approach (for details, see [51]).

We analyzed the correlations of the three connectivity metrics for three different scenarios to characterize movements with strategic significance in terms of intra- and inter-subpopulation movements. The three scenarios include movements within the eastern (EE) and western (WW) subpopulations and between subpopulation edges (WE)—i.e., border of each subpopulation closest to the other. We considered only the edges and not subpopulation centers to make a fair comparison between methodologies since Linkage Mapper and Conefor process the links of spatially adjacent nodes and only considered movements between subpopulation edges. Another reason is that movements inside subpopulations are already considered in two scenarios and happen in suitable habitat areas where these movements are feasible and frequent. In addition, connectivity conservation efforts are mostly focused on the low permeability gap between subpopulations [38,51].

3. Results

3.1. Habitat Suitability Models

3.1.1. Operational Scales and Predictors Effect

The operational scales of the vegetation variables selected by the univariate models varied among models. Foraging resources were similar between global (0.5 km) and continental (2 km) datasets but differed significantly from the national one (16 km), which was much broader. Forest area, on the contrary, was analogous for the three datasets—variation between 0.5 and 1 km. Terrain ruggedness influenced at a relatively small scale of 2 km and all human pressure variables best performed at 16 km in the three models (

Table 2. Habitat suitability models for the three datasets resulting from different data sources used to calculate vegetation variables (i.e., foraging resources and forest area). Note that human pressure and terrain ruggedness variables are constant in the three models. “Scale” is the operational scale (km) for each environmental variable and data source. Significance levels: * <0.001. See Figure A1, Figure A2 and Figure A3 for the layers of the operational scales of the variables.

Predictor	Global		Continental		National
	Scale	Coefficient	Scale	Coefficient	Scale	Coefficient
For.Resources	0.5	1.081 *	2	0.547 *	16	1.854 *
Forest Area	1	−0.512 *	1	0.210 *	0.5	0.400 *
Ruggedness	2	1.048 *	2	1.032 *	2	0.982 *
Buildings	16	−1.324 *	16	−1.775 *	16	−1.726 *
Highways	16	−1.376 *	16	−1.315 *	16	−0.156
Roads	16	−2.004 *	16	−1.580 *	16	−0.523 *
AUC (95% CI)	0.900 (0.896–0.904)		0.908 (0.904–0.911)		0.936 (0.933–0.939)

).

Vegetation variables had a positive effect (Table 2) on habitat suitability in most of the cases, especially the foraging resources of the national dataset. Together with foraging resources, terrain ruggedness was the environmental variable that most fostered bears’ habitat suitability in every multivariate model. Forest area had a negative influence on the global model, while it was slightly positive in the continental and national models. In all cases, the relation of human pressure variables was unfavorable.

3.1.2. Model Performance and Habitat Suitability Spatial Pattern

All the models yielded a strong predictive performance with high AUC scores (Table 2). The national model achieved the highest AUC (0.936) while the lowest belonged to the global model (0.900). Results of DeLong’s test showed a significant difference between the three models (p-value < 0.05). Overall, the spatial pattern of the areas with high habitat suitability was very similar and almost identical between the global and the continental models, although with differences in their sharpness. High values appeared in two cores, which matches the brown bear distribution range in the study area. Both cores were separated by low-suitability areas where two main highways occur (AP-66 and N-630; see Figure 2). High suitability also appeared towards the southwestern and southeastern borders of the study area in the global and continental model predictions.

3.2. Connectivity Analysis

3.2.1. Corridor Comparison

Results from the connectivity analysis were variable according to the methodology used to compare the different datasets. Links obtained through least-cost modeling showed mostly similarities between datasets and the genetic model, especially inside both subpopulations (see Figure 3 and Figure A4 in Appendix B for the genetic model). For the connections between subpopulations, the links derived from three datasets showed few differences as they mostly concurred, and when they did not, their routes were often parallel. Nevertheless, their differences were more evident for longer links in the east border of the eastern subpopulation and in the northwest.

The three datasets revealed moderately different current flow patterns when using circuit theory (see Figure 4). The global model drew a sharper and more scattered pattern of corridors, in accordance with its corresponding habitat suitability map. On the contrary, the spatial pattern of the continental model was more homogeneous and concentrated in fewer corridor areas. The national model showed a less dendritic pattern than the global model but with less condensed flow values than the continental model. In the three models, high current flow values appeared in similar areas but with more irregular patterns within the eastern subpopulation. Current flow values were low between subpopulations in the three models.

3.2.2. Comparison of Connectivity Metrics

We found high correlations when comparing connectivity metrics among the three resistance surfaces (see

Table 3. Pearson’s correlation of connectivity metrics: effective (least-cost modeling) and resistance (circuit theory) distances, and difference in probability of connectivity index (dPC) values (PC index) between the paired combinations of global, continental and national models and of each of these three models against the genetic model for each subpopulation (subp.) and between subpopulation edges.

Pearson Corr.	Within West subp.	Within East subp.	Between subp. Edges
Least-cost modeling
Global–National	0.985	0.953	0.986
Global–Continental	0.996	0.988	0.995
Continental–National	0.987	0.965	0.995
Global–Genetic	0.951	0.919	0.972
Continental–Genetic	0.966	0.906	0.985
National–Genetic	0.972	0.920	0.994
Circuit theory
Global–National	0.897	0.765	0.735
Global–Continental	0.961	0.930	0.939
Continental–National	0.905	0.842	0.750
Global–Genetic	0.747	0.649	0.582
Continental–Genetic	0.742	0.640	0.494
National–Genetic	0.748	0.663	0.774
PC index
Global–National	0.214	0.610	–0.094
Global–Continental	0.427	0.656	–0.044
Continental–National	0.690	0.640	0.490
Global–Genetic	0.313	0.097	0.059
Continental–Genetic	0.218	0.121	0.087
National–Genetic	0.410	0.394	0.678

). Under least-cost modeling, Pearson’s correlation estimates (r) of effective distances were remarkably strong, with coefficients higher than r > 0.95 between global, continental and national models. The highest correlation recorded was inside the western subpopulation between the global and continental models (r = 0.996). Strong correlations remained when comparing the effective distances of the three models with the genetic model (r > 0.9).

In comparison to least-cost modeling, the correlation patterns were similar when using circuit theory, although weaker correlations were observed. The best correlations remained between the global and continental models for every type of movement. The correlation coefficients were moderate (0.8 < r < 0.5) between each of the three models and the genetic connectivity model. The highest coefficients were recorded by the national model, especially for movements between subpopulation edges.

The dPC values showed less consistent correlations among models (see Table 3). Continental and national models showed moderate correlations. However, the correlations of the global model were weaker, with poor correlations between subpopulation edges. The global and continental models barely correlated with the genetic model. Conversely, the national and the genetic dPC values correlated moderately for movements between subpopulations.

4. Discussion

The choice of spatial data source is an important step in ecological modeling since the resolution of data is a key factor for landscape definition that deeply influences the nature of the results [60,61]. Accordingly, our results show similar outcomes but also meaningful divergences when modeling habitat suitability and connectivity using different data sources.

4.1. Habitat Suitability Models

Our results show similar predictive performances for the three habitat suitability models, which suggest that Sentinel-1 and 2 imagery and CLMS pan-European products are solid alternatives for modeling habitat suitability. Besides minor differences in AUC scores, habitat suitability maps showed analogous spatial patterns of high-suitable areas with the three models correctly predicting the distribution range of the bear in the Cantabrian Range.

The fine-spatial resolution data included in every model—i.e., Sentinel-1 and 2, FTY and LiDAR—may be the main cause behind the good performances, since finer grains make land cover classes more distinguishable for the model [62,63]. The high spatial resolution led to the selection of small operational scales (0.5–1 km) in the three models for the forest area variable and allowed the identification of isolated high-suitability areas (forest patches with abundance of food resources) in such a heterogeneous landscape [16]. However, the delineation of high-suitability areas did not have a detailed pattern in the continental and national maps. This could be due to CLC18 and FMS as their broader scale (MMU of 25 and 2.25 ha, respectively) may have affected the detail of the results.

To match the resolution of the bear presence data, we resampled the resolution of the environmental variables to 1 ha, coarsening the fine pixel size of some layers (e.g., from LiDAR, FTY and SLCM). Brown bear locations had to be resampled to 1 ha to avoid pseudoreplication since multiple traces could have been caused by the same individual. GPS and tracking data could improve the accuracy of the locations, allowing using the abovementioned environmental variables with their best spatial resolution. This could give new insights into the influence of high spatial resolution on habitat suitability modeling. However, for most cases, brown bear habitat preferences are more influenced by coarser operational scales and not by the finest scale available of 250 m (Table 2). Despite the operational scales, the high spatial resolution of the layers did influence how the configuration of the landscape was characterized, considering small habitat patches and landscape elements that would have been ignored with coarser spatial resolutions before resampling.

Thematic resolution may explain the slightly different predictive performances among models, as performance increased when their datasets classified more vegetation types [12,62]. For instance, the identification of multiple vegetation types included in the national model (due to the FMS) allowed characterizing important foraging resources for the brown bear—in comparison to the SLCM which only classified one class of forest. The richer thematic resolution of the national dataset may also explain the different operational scales selected for the foraging resources variable.

The inaccuracies of input layers could also produce lower performances in the continental and global models. Both the FTY (minimum overall accuracy of 0.90) [25] and the SLCM layers (kappa index of 0.89) [27] had high accuracies. However, some misclassifications are also visible, especially in the SCLM. The steep slopes negatively affected Sentinel-1 data, producing classification errors in mountainous areas where vital natural resources for the brown bear are abundant. More misclassifications appeared between forested areas and shrublands, both important for the definition of foraging resources [27]. These inaccuracies could explain the negative influence of the forest area variable in the global model (Table 2) as some relevant suitable areas for the bear were not correctly represented as forest in this model. Sentinel-1 allows mapping vegetation variables regardless of cloud conditions that limit optical imagery [64], but the potential issues derived from its use in hilly landscapes should be addressed [65].

Our findings support the capabilities of satellite imagery to predict species habitat preferences in heterogeneous landscapes [13,66]. However, this may not be appropriate for all circumstances. The considered global and continental datasets may be more suitable for cases with broad ecological requirements, as these may not attain enough detail in identifying a wide variety of natural resources, which can be important for understanding species’ habitat requirements at small scales [32,67,68]. When trying to understand specific habitat–species relations, spatial analysts should use data that enrich the thematic resolution [69] and that complement the high spatial resolution of the input vegetation variables. For this purpose, free remotely sensed data sources, such as the recent space-borne NASA LiDAR, can be useful to inform about the structure of vegetation [22].

4.2. Connectivity Models

The similar trajectories, current flow patterns and high correlations of distances among models when applying least-cost and circuit theory methodologies indicate that varying data sources did not profoundly affect connectivity results. However, even the slight variations recorded in the results could have considerable implications for conservation [51,63]. For instance, lower correlations in resistance distances and slightly different corridor patterns suggest that the data source choice can have implications for corridor identification through circuit theory. Much of the difference among current maps may be derived from the effect of low-spatial resolution data (CLC18, FMS; see Section 4.1). Furthermore, the moderate to low correlations of dPC values (Table 3) indicate different contributions to connectivity for analogous links. Thus, each connectivity model would prioritize distinct corridors, which could have repercussions in the landscape planning of connectivity measures. In any case, predicted corridors at regional scales should be used as a first step in characterizing connectivity networks and in combination with studies conducted at different scales. More detailed data suitable for local scales should be used to ascertain if the predicted corridors are functional and better define the routes, habitat patches and landscape features before implementing specific measures of conservation projects [70,71].

Grain size can slightly but importantly affect the total cost of movement of a corridor [50,72], as finer grains capture landscape and vegetation patterns that should be represented to not misestimate corridor trajectories [73]. Considering the similar spatial grain of the vegetation variables, the differences among connectivity model outcomes may be due to variable thematic resolutions for estimating their resistance surfaces [72].

Genetic connectivity models are more representative of dispersal movements and successful reproduction between subpopulations, making them more appropriate for corridor planning [50,51,74]. Thus, the strong correlations (Table 3) with the genetic model when using least-cost modeling suggest that the spatial datasets considered in this study can offer robust connectivity results under this approach. It should be noted that least-cost modeling is deterministic and predicts only optimal paths [51,55,75]. Despite having similar effective distances for analogous paths, minor changes in their trajectories could focus the implementation of conservation resources on different places. This effect was more prominent as the Euclidean distance between nodes increased (Figure 3). Under circuit theory, the correlations with the genetic model were moderate but the strong coefficient between subpopulation edges of the national model highlights this model as the best option to predict corridors through degraded landscapes. A similar result was obtained for dPC values, supporting the national model also as the best alternative to assess the corridors’ importance between habitat areas. These results are particularly relevant, considering that conservation efforts often focus on the identification of corridors between disconnected habitat patches to enhance gene flow, species dispersal and adaptation to current threats [51].

4.3. Other Considerations on the Choice of Spatial Data Source

Above, we discussed the effects that the choice of spatial data source can have on two types of ecological models (i.e., habitat suitability and connectivity). If multiple data sources are available, our results advocate for the data source with the highest spatial and thematic resolution. However, the regional availability of data can limit the choice of data source, which in most cases would be restricted to satellite data, such as Sentinel-1 and 2 imagery. Nevertheless, the temporal resolution of data can also be essential, especially for dynamic habitat and connectivity studies [76]. Satellite remote sensing vegetation variables offer a great advantage as these can be updated regularly with a high frequency—5 to 6 days of average revisit time intervals for Sentinel sensors at the equator [20,21]. Despite being derived from satellite data, CLMS products do not offer frequent updates—every 6 years for Corine Land Cover and every 3 years for the FTY. This period can be even longer for nationally restricted data; e.g., Spanish LiDAR data are updated approximately every 7 years, while the FMS has not yet been updated. Consequently, the temporal resolution would favor the choice of satellite imagery to model ecological processes when the time factor is important. Additionally, the combination of satellite data with other data sources can be recommendable to offset the caveats that these sensors can have [22]. Nonetheless, satellite imagery processing requires great expertise and effort to map environmental variables. This limitation supports the choice of already processed cartography, such as CLCS products, which allow conservationists to focus almost directly on the modeling process [77].

5. Conclusions

Our results support free satellite imagery as a global and powerful source to estimate vegetation variables for habitat suitability and connectivity modeling. Already processed land cover products (such as the ones from CLMS) are also an effective alternative for studies focusing on the European continent. Therefore, widely available data sources are suitable to produce meaningful results that guide the design of consequent conservation measures and nature management plans. However, some considerations depending on the required spatial, temporal and thematic resolution need to be taken into account. Globally available satellite data have the advantage of being updated frequently, but they may require sufficient processing skill and can be time-demanding and computationally intensive. Satellite data may be limited in many cases to understand precise habitat–species associations, which requires high thematic resolution besides fine spatial detail. Failing to map a sufficient diversity of natural resources for a species can have detrimental effects in some cases, such as in the identification and prioritization of corridors. For this reason, more effort should be applied to improve the thematic resolution of vegetation variables coming from broadly accessible sources, which can be achieved by combining sensors and other open source spatial data sources. The comparison presented by this study offers valuable insights to better understand the potential implications, biases and limitations that spatial data source selection has on the results of habitat suitability and connectivity modeling, which is key for the correct definition of conservation measures.