Biodiversity and infrastructure interact to drive tourism to and within Costa Rica.

SignificanceTourism accounts for roughly 10% of global gross domestic product, with nature-based tourism its fastest-growing sector in the past 10 years. Nature-based tourism can theoretically contribute to local and sustainable development by creating attractive livelihoods that support biodiversity conservation, but whether tourists prefer to visit more biodiverse destinations is poorly understood. We examine this question in Costa Rica and find that more biodiverse places tend indeed to attract more tourists, especially where there is infrastructure that makes these places more accessible. Safeguarding terrestrial biodiversity is critical to preserving the substantial economic benefits that countries derive from tourism. Investments in both biodiversity conservation and infrastructure are needed to allow biodiverse countries to rely on tourism for their sustainable development.

The tourism sector is well-poised to generate win–win approaches to biodiversity conservation and sustainable development given its nonextractive nature and its dependence on scenic beauty (1). However, for sustainable tourism to succeed as a strategy for biodiversity conservation, the role that biodiversity plays in driving tourism patterns needs to be better understood. On one hand, wildlife and nature motivate a significant portion of global tourism (2), and protected areas with higher species richness tend to attract more tourists and yield higher economic benefits (3). On the other hand, tourism hotspots also tend to occur in places where more human-built infrastructure (e.g., hotels, roads, and airports) enables access (4, 5). Studies have reached mixed conclusions on the relative importance of biodiversity and accessibility for tourism, and little is known about how they work in concert (6–11). Given the potential negative impacts of infrastructure on biodiversity conservation, their relative contributions to tourism deserves explicit study, particularly in developing countries where both biodiversity conservation and sustainable development are urgently needed (7, 12).

Drivers of tourism patterns across landscapes have been explored through questionnaire surveys and structured interviews that ask tourists about their affinity for landscape features (13–15) and through spatial models that predict recreation using photographs (e.g., geographically weighted regression, MaxEnt) (16, 17). Recently, geo-tagged photographs and species lists shared on social media platforms have become popular tools for tourism-focused research (18–22). These studies, however, typically focus narrowly on the role of single taxa (22, 23) and landscape attributes without accounting for species diversity (14, 20), or they focus only on the role of infrastructure as a driver of tourism (5). Recent advances in satellite Earth observations make it possible to capture more of the ecosystem heterogeneity that can drive variability in species distributions, compared to more conventional modeling based on land cover (24–26). An integrated approach is needed, linking species richness of multiple taxa along with infrastructure variables, both modeled and mapped through high resolution Earth observations. Such an approach could be scaled up to larger regions and applied globally, helping to identify where biodiversity is playing a significant role in driving tourism, such that governments and the tourism sector can prioritize investments in biodiversity conservation.

Here, we ask three questions in the iconic case of Costa Rica. First, what is the contribution of species richness (of vertebrate taxa) in explaining patterns of tourism in protected areas and also country-wide? Second, how similar are the patterns for birdwatching tourism compared to those of overall tourism? Third, where in the country is biodiversity contributing more than other factors to birdwatching tourism and to overall tourism? We predict that vertebrate species richness is more important for driving tourism in protected areas than in the rest of the country, because nature-seeking tourists often go to protected areas to find wildlife (3). We also expect that birdwatching tourism is predicted by richness of threatened and endemic bird species rather than total species richness, given birdwatchers’ preferences for rare birds (27). We predict a saturating relationship between species richness and tourism, because beyond a large number of species additional species are unlikely to contribute more to tourism (28). Finally, we predict that national-level tourism is better explained by infrastructure (such as roads and hotels) and distance to water than by biodiversity, because tourists going to Costa Rica often seek activities such as surfing and relaxing in beach resorts (29, 30). We predict nonlinear effects of proximity to roads and water, because a place is deemed inaccessible if it is further away from roads, and a beach tourist destination is also either close to water or not a destination at all. Access diminishes rapidly over a few miles (31).

To answer our research questions, we analyze patterns of tourism at two different spatial scales. First, we analyze tourism patterns in protected areas only. Many protected areas provide reliable data on visitation rates and biodiversity tends to be higher inside than outside protected areas (29). However, given that they are often visited by tourists who are already interested in biodiversity, protected areas are not representative of all tourism patterns (30). Second, we evaluate the relative importance of biodiversity to all tourism across Costa Rica (excluding offshore islands). Investigating tourism across the whole nation may give a better understanding of how biodiversity contributes to tourism writ large and not only for tourists with a predisposition for finding wildlife. At both the protected area and national scales, we use a modified MaxEnt model that integrates species distribution models for 66 terrestrial vertebrate species (including amphibians, reptiles, and mammals) and 699 bird species, based on remotely sensed climate and habitat variables, with spatial patterns of infrastructure (hotel density and distance to roads) and distance to water. We measure tourism in two ways for both scales: using eBird checklists as a proxy for birdwatching (32), and using Flickr photographs as a proxy for all international and domestic tourism (19).

Costa Rica is an ideal country to explore these questions because tourism represents 7% of the national gross domestic product and employs 3% of the working population directly and a further 9% indirectly (33). Approximately 70% of all international visitors to Costa Rica state that the wildlife, dramatic scenery, and opportunities for adventure sports are the main motivation for visiting the country (33, 34). However, the importance of biodiversity as a factor that influences tourism to, or domestically within, the country has not been evaluated (apart from very local studies) (35, 36). The Central Bank of Costa Rica is currently piloting a nature-based tourism account under the United Nations System of Environmental Economic Accounts (UN SEEA). The state of the art with this methodology is to attribute value to different ecosystems, which may vary widely in their biodiversity. Understanding the relationships between biodiversity and tourism in Costa Rica is a key step toward maintaining the vibrant ecotourism industry and can serve as an example for other biodiverse nations that often look to Costa Rica as a leader in sustainable development (37).

Results

Tourism inside Protected Areas and Country-wide.

As expected, vertebrate biodiversity (species richness of amphibians, reptiles, and mammals) matters more to tourism within protected areas than country-wide. The spatially explicit model aimed at explaining tourism within protected areas, as measured by density of Flickr photographs-user-days (PUDs) per protected area shows a better fit (area under the curve [AUC] = 0.680; —Model 5) than the nation-wide model (AUC = 0.608; —Model 10). The protected area model shows that protected area size (50.8%), hotel density (27.1%), distance to roads (10.3%), and distance to water (11.1%) are the most important predictors in terms of permutation importance (—Model 5). We find negative effects of protected area size, distance to roads, and distance to water, suggesting that larger parks and parks further from roads and water have lower densities of tourists compared to smaller parks with greater accessibility (Fig. 1 and ).

Fig. 1.

Biodiversity has a positive effect on tourism inside protected areas. Bird species richness (A), richness of threatened and endemic bird species (B), and richness of amphibians, reptiles, and mammals (C) all have a nonlinear effect on tourism. Birdwatching and overall tourism decline when distance to roads increases (D), and with increased protected area size (E), as well as increased distance to water (F). Protected areas with intermediate levels of hotel density experience the highest tourism (G). These curves show the results of ten k-fold replicated MaxEnt runs using either Flickr PUDs (solid lines) or eBird CUDs (dotted lines) as the response variable and the corresponding environmental variable as the predictor. The line represents the mean of the ten replicates, while the shaded region is drawn between the minimum and maximum prediction of the 10 replicates.

The species richness of amphibians, reptiles, and mammals, is the least important predictor with the lowest values of permutation importance (0.6%) and shows a curved pattern, suggesting that parks with intermediate levels of vertebrate species richness have more tourists, though effects are muted (). Importantly, the role of amphibians, reptiles and mammals is larger in explaining tourism inside protected areas compared to just bird species richness (0.4% permutation importance, —Model 3, AUC = 0.679), and much lower compared to the role of richness of threatened and endemic birds (10%, —Model 4, AUC = 0.690). Nonetheless, in the three models, the effect of biodiversity indicates that tourism is higher with intermediate levels of diversity (). In all of the models, the main predictor is protected area size with a negative relationship (—Models 3–5).

Analyses at the country-wide scale reveal that accessibility is the main driver of tourism. Indeed, the main predictors explaining tourism patterns, as measured by Flickr PUDs, are distance to roads (∼50% of permutation importance in all three models, —Models 8–10), and hotel density (∼20% of permutation importance in all models, —Models 8–10). Distance to roads also has a negative effect, reaffirming that places closer to roads have more tourists. Hotel density has a quadratic effect, suggesting that people visit areas with intermediate levels of hotel density (Fig. 2 and ). Distance to protected areas is the third most important predictor across models (∼20% of permutation importance in all three models) and it exhibits a negative linear effect, indicating that being close to a protected area yields more tourism.

Fig. 2.

Across Costa Rica we find that tourism increases when species richness of birds (A), richness of threatened and endemic birds (B), and richness of amphibians, reptiles, and mammals increase (C). Tourism is predicted to decrease as distance to roads (D), distance to protected areas (E), and distance to water increases (F). The relationship between hotel density and tourism shows a nonlinear pattern, where intermediate levels of hotel density lead to highest levels of tourism (G). These curves show the results of ten k-fold replicated MaxEnt runs using either Flickr photographs (solid lines) or eBird CUDs (dotted lines) as the response variable and the corresponding environmental variable as the predictor. The line represents the mean of the ten replicates, while the shaded region is drawn between the minimum and maximum prediction of the ten replicates.

In keeping with the protected area analysis, the biodiversity variable that better predicts overall tourism at the country-wide scale is the richness of threatened and endemic birds (5.3%; AUC = 0.612, —Model 9). The importance of amphibians, reptiles, and mammals (0.1%, —Model 10) and bird species richness (3.8%, —Model 8) are lower. Nonetheless, in all models, we find a positive effect of biodiversity in explaining tourism (Fig. 2 and ). Last, distance to water has a permutation importance ranging from 1.9 to 3.1% and has a negative nonlinear effect, suggesting that people prefer being closer to water, within a few kilometers and generally not more.

Birdwatching Tourism inside Protected Areas and Country-wide.

Birdwatching tourism is patterned similarly to overall tourism but exhibits greater nonlinearities with regards to infrastructure and biodiversity variables. Models explaining birdwatching inside protected areas are the best performing models across our entire analysis, across the two different scales and two types of tourism (AUC = 0.776 for the model with all birds, and AUC = 0.780 for the model with threatened and endemic birds, —Models 1, 2). The main predictors for birdwatching inside protected areas are distance to roads (61.1% and 54.2% for the model with bird species richness [Model 1] and for the model with richness of threatened and endemic birds [Model 2] respectively), and protected area size (13.1% and 14.8%), both showing negative nonlinear effects (Fig. 1 and ). We find that the richness of threatened and endemic birds (7.1%) is a better predictor of birdwatching compared to all birds (2.1%). We find that birdwatching changes nonlinearly with respect to bird species richness, showing that birdwatching tourism peaks at ~200 species and at ~100 species of threatened and endemic birds (Fig. 1).

Our model examining birdwatching tourism at the country level, as measured by eBird checklists with bird richness as the biodiversity predictor variable, has a similar fit compared to the model with threatened and endemic bird richness (AUC = 0.723, Model 6 vs. AUC = 0.721, Model 7). Distance to roads (∼58%) and hotel density (∼13%) are the main predictors of birdwatching tourism nationally (—Models 6 and 7). We find stronger negative relationships between predictor variables and birdwatching tourism than overall tourism (Fig. 2). Specifically, we find negative quadratic effects for distance to roads () and higher negative slopes of distance to water () compared to birdwatching at the protected area level (), suggesting that being close to a road or to water leads to more birdwatching up to certain distances (<50 m for roads and <50 m for water) before dropping significantly (Fig. 2). Conversely, we find a positive saturating effect for hotel density, indicating that places with more hotels have more birdwatching but saturate at 0.75% of hotel density. Species richness of threatened and endemic birds is the third most important predictor in explaining birdwatching tourism at the country-wide scale (12.1%, —Model 7). Importantly, we find that the richness of threatened and endemic species is a better predictor of birdwatching patterns compared to all birds (10.3%, —Model 6).

Contribution of Biodiversity to Spatial Patterns of Tourism.

Applying our spatial regressions to national maps clearly demonstrates the importance of Costa Rica’s major population centers and road networks to tourism patterns (areas in dark green in Fig. 3 ). When the models are run with and without biodiversity, new patterns emerge. Patterns are stronger for birdwatching tourism (Fig. 3) than for overall tourism (Fig. 3), but in both cases, without considering biodiversity, the tourism value of the Cordillera de Talamanca (dark red areas in Fig. 3 ) is underestimated, while the tourism value of the Nicoya Peninsula is overestimated (blue areas in Fig. 3 ).

Fig. 3.

Tourism is higher in places with both biodiversity and infrastructure. Maps indicating higher (green) and lower (white) tourism prevalence standardized between 0 and 1 (A, D). Maps showing the predicted birdwatching tourism at the country-wide scale as predicted by eBird checklists (A) and overall tourism as predicted by Flickr photographs (D). Difference maps showing how excluding either richness of threatened and endemic birds affects tourism as measured by eBird CUDs (B) and Flickr PUDs (E) or how excluding infrastructure (hotel density and distance to roads) affects tourism (C based on eBird CUDs, F based on Flickr PUDs). Dark red areas (B, C, E, F) indicating that tourism is underestimated if those variables are excluded vs. blue areas indicating that tourism is overestimated if those variables are excluded. The dotted lines in all panels represent the protected area network, with the largest park (La Amistad International Park) shown in the south-east of the country.

Models with and without infrastructure show that tourism in intact tracts of forest is overestimated, such as the Amistad International Park (dark blue area in the southeast), the Corcovado National Park in the Osa Peninsula (protected area in the southwest), and the Barra del Colorado Wildlife Refuge (northeast of the country), whereas areas along the Pacific coastline and near urban centers tourism are underestimated (Fig. 3 ).

Results show that tourism patterns are higher when both biodiversity and infrastructure are considered together, especially in areas near Monteverde, and the Volcán Poás National Park (Fig. 3, center of the country). We find some spatial mismatches between tourism and birdwatching; in particular, birdwatching is more reliant on the road network than overall tourism, and certain protected areas are more important for tourism than for birdwatching, such as the Amistad International Park and the Barra del Colorado Wildlife Refuge (Fig. 3 ).

Discussion

Our results indicate that if Costa Rica invests in both infrastructure development and biodiversity conservation, they will continue deriving economic benefits from their natural assets (ceteris paribus). We show that while the role of infrastructure, such as hotels and roads, is larger in driving tourism compared to the role of biodiversity, tourism is highest in places where both biodiversity and infrastructure are present, particularly in mountainous areas. It is well-known that the tourism footprint is strongly associated with access and is particularly dependent on the presence of roads (31, 38, 39). Our findings reaffirm that proximity to roads increases tourism inside protected areas and country-wide. Even for birdwatchers, being close to roads is a better predictor of where they go rather than the richness of threatened and endemic birds. Thus, our results echo prior studies, finding that accessibility is more important than the ecological value of a place in driving tourism, including nature-based tourism (40). However, failing to adequately account for nature’s contributions will overestimate the value of highly accessible but ecologically depauperate places, and will underestimate the importance of sustainable development to maintaining current levels of benefits.

Distance to water is less important in driving tourism patterns than we expected. Considering that ∼70% of all international visitors to Costa Rica list ecotourism, adventure, and beaches as their main motivation for visiting the country (34), we expected proximity to water to play a larger role in driving tourism country-wide than our analyses suggest. While beach-goers and surfers constitute a significant portion of tourism, many visitors are seeking other activities such as birdwatching that are not dependent on proximity to water. These results might indicate that tourist activities in Costa Rica are more diverse than previously reported (41, 42).

Aligned with our predictions, we find that the richness of amphibians, reptiles, and mammals explains tourism inside protected areas better than country-wide. However, the contribution of this variable to tourism is smaller than we expected. The slightly bigger predictive power of these vertebrates, compared to birds, of tourism inside protected areas might be explained by visitors’ preferences toward charismatic fauna, such as monkeys (e.g., White-faced capuchins, howlers) and spectacled caiman (Caiman crocodilus), which are some of Costa Rica’s flagship species and main attractors of tourists (43, 44).

Moreover, while the slope of the relationship is fairly shallow, we see that intermediate levels of diversity yield the highest tourism; this echoes results from Namibian conservancies, where more diverse ecological communities provide tourists with a greater range of viewing opportunities (45). With respect to birdwatching, it is unsurprising that biodiversity is more important in explaining this form of nature-based tourism, compared to all tourism (22). It is important to recognize that the nonlinear effect of species richness on birdwatching and overall tourism may indicate that the identity of the birds in communities is more important than total richness. These results are consistent with studies showing that species richness alone does not drive cultural value, as in the case of wildlife viewing and preferences for beautiful landscapes with flowering plants (11, 15). Many birders seek out rare and endemic species, whose protection may depend on greater conservation effort that, in turn, may bring greater economic benefit.

Our spatial models identify some differences between where birdwatchers go and where general tourists go. While spatial congruence of species richness tends to be high between vertebrate taxa (e.g., places with more birds also have more mammals), it is not always congruent with the richness of rare and endemic vertebrate species (46). Thus, in our study we show that the highlands of Costa Rica have more birdwatching tourism than any other place. Indeed, the highlands are where the endemic birds are more common (47), and where the most preferred species by birdwatchers can be found (e.g., Resplendent Quetzal and Ornate Hawk-Eagle) (48) (Fig. 3).

Moreover, we find that while birdwatching, like all tourism, increases with infrastructure, it is highest where both biodiversity and infrastructure co-occur. Thus, if biodiversity is not considered, places like Amistad International Park, and the Corcovado National Park in the Osa Peninsula would be predicted to have less tourism than currently observed (Fig. 3). We conclude that a degradation of the country’s ecosystems could negatively impact the country’s economy. In typical years, Costa Rican tourism employs 160,000 people directly and 450,000 more people indirectly (33)—although we recognize that the COVID-19 pandemic has severely impacted the tourism sector and the local communities that rely on tourism in the past years (49). Loss of biodiversity could mean further attrition in tourist’s visits, with cascading impacts to revenues and jobs.

Funding from multilateral development banks, bilateral cooperation agencies, and the private sector could prioritize investing in infrastructure to support sustainable tourism in places where biodiversity is high (50) and to advance the international policy agendas of mainstreaming biodiversity in the tourism and infrastructure sectors (51, 52). Nonetheless, it is critical to recognize that increased access to protected areas by expanding road networks in undisturbed areas (e.g., Amistad International Park), might cause biodiversity declines, as has been demonstrated in other tropical regions (53, 54). Therefore, a crucial management challenge for Costa Rica, and other countries depending on tourism, is to plan infrastructure development in the least destructive ways. This means not only avoiding the most sensitive areas and maintaining the balance between an ecosystem’s integrity and its accessibility to the public (55), but also policies, plans, investments, and cultural norms reinforcing the conservation of nature (56).

The relationship between tourism, biodiversity, and infrastructure documented here may be sufficiently general to translate to other developing countries with similar species richness, but the relative importance of biodiversity and infrastructure for tourism might vary across spatial scales. For instance, for nature-based activities such as birdwatching and wildlife viewing, the presence of certain species and the diversity of species might determine where people go (15). Birdwatchers, in particular, are known to be collectors who prefer endemic and range-restricted species (57), which might factor into the decision of which country to visit, whereas infrastructure such as hotels and roads might explain where people go birdwatching once they are in a country (5).

We also note that Costa Rica is politically stable and has elected a series of governments that prioritize environmental management and protection (58). Governance indeed is a key factor determining tourism success, even inside protected areas (59). Thus, tourism in other places that have high biodiversity and high road density might not yield the same economic revenues due to the lack of political stability and security. Moreover, leakage can occur when some areas are developed for tourism while others with similar characteristics get degraded due to the lack of oversight and regulation (60). Future studies could compare the role that biodiversity and infrastructure play in driving tourism relative to the role of governance in locations where tourism is a promising solution for sustainable development, but political conditions are less stable (e.g., Colombia and the Democratic Republic of Congo) (61, 62).

Our research builds on a growing number of studies showing that social media are an effective source of information on tourism (19, 63, 64). While MaxEnt was originally created to model the distributions of species based on environmental predictors such as climate and land use (65), it has recently been used for mapping social values for ecosystem services (16, 17). To date, MaxEnt applications involve either the mapping of cultural ecosystem services by identifying people’s preferences for landscape attributes (16), or species distribution modeling for conservation planning (66). They are rarely done together, except in very small-scale studies (67). Therefore, the methodological advancements presented in this paper link Earth observations (via land-use maps derived from satellite imagery) with vertebrate species distributions, and with social media data at national and subnational scales.

While this method is powerful for informing national planning decisions, particularly by highlighting the places with high tourism demand and the places where biodiversity has a disproportionate effect on tourism, it also has some limitations. One is the potential taxonomic bias in the Global Biodiversity Information Facility (GBIF) dataset used to generate our species distribution models. Specifically, many rare species may not have sufficient observations in GBIF to meet our stated requirements for inclusion (68). For example, we were able to model many more birds than amphibians, reptiles, and mammals, and yet tourist’s surveys indicate that these latter three taxa are their preferred animals to find (43). Therefore, our model is likely underestimating the role of rare species, such as jaguars, in driving tourism (69). Future studies are needed to evaluate the role of rare species on driving tourism, perhaps with complementary research methods (e.g., choice experiments).

In conclusion, by using Costa Rica as an example, we provide a pathway for future evaluations of the role of biodiversity in driving tourism and nature-based tourism at national scales, which is the most relevant scale for development planning and biodiversity conservation strategies (70). Our proposed methodology and approach could inform UN SEEA accounts and other approaches aimed at mainstreaming biodiversity in the tourism sector and scaling it up to other countries (71).

Materials and Methods

Study Area.

Costa Rica’s land surface covers 51,100 km2 and harbors ∼5% of the world’s macroscopic species, making it the world’s number one country in terms of species density per 1,000 km2 (72). This small country encompasses 95 distinct climatic zones (73), ranging from the submoist dry and very warm climate of the Nicoya Peninsula lowlands in the northwest to the very wet and cold climates in the Cordillera de Talamanca, which runs from the northwest to the southeast of the country and rises to an elevation of 3,819 m (73). Costa Rican climate is influenced by both the Atlantic and the Pacific oceans. Its mean annual temperature varies from 26.7 °C on the north Pacific coast, to 26 °C on the Caribbean coast, to 6 °C on the highest peak of the mountain range. Total annual precipitation varies from 1,300 mm to 7,467 mm, and average monthly relative humidity varies between 65 to 90% (73). The climatic variation results in a tremendous diversity of ecosystems. The most recent classification suggests 14 major ecosystem types that range from the páramos (above treeline), to seasonal tropical dry forests, evergreen tropical rainforests, mangroves, rivers and estuaries, bogs, marshes, and swamps (74).

Data Collection.

To test the relationship between predictor variables and tourism and birdwatching we conducted analyses at two spatial scales. First, we evaluated such relationships at the country level, and second, we did so at the protected area level (n = 108 protected areas, ). We collected data for eight predictor variables (1-Distance to roads, 2-Hotel density, 3-Distance to water, 4-Distance to protected area, 5-Protected area size, 6-Richness of amphibians, reptiles and mammals, 7-Bird species richness, 8-Richness of threatened and endemic birds) and two response variables (1-Flickr photographs, 2-eBird checklists) (). What follows presents the methods for calculating each predictor variable.

Distance to roads.

We extracted the roads in Costa Rica from the Global Roads Open Access Database in 2018. This layer was then rasterized at a 915 m resolution (in Universal Transverse Mercator [UTM] zone 16N), with all pixels intersecting roads being assigned a value of 1. We then calculated the distance to roads using the distance command in the R package raster. We used rasters of 915 × 915 m of resolution throughout the study.

Hotel density.

We first downloaded a list of hotels registered in Costa Rica by 2018 through the GeoNames database (https://www.geonames.org/). This list contains the latitude and longitude of hotels, nature lodges, and hostels throughout the country (). We converted this layer to UTM zone 16N coordinates. We used the “kde2d” function from the MASS package in R (49) to perform a kernel density calculation. Using a Gaussian decay function with a SD of 25 km we generated a hotel density map for the entirety of Costa Rica matching the resolution of our other rasters (915 × 915 m).

Distance to water.

We created a composite layer of water bodies in Costa Rica by combining shapefiles of inland rivers and lakes with a global ocean shapefile, both downloaded from the Natural Earth database (https://www.naturalearthdata.com/). We converted this shapefile to a raster matching the resolution of the previous layers. All pixels that touched polyline features (rivers), and pixels that were covered by at least 50% of a polygon (lakes and oceans) were converted to values of 1 (indicating presence of water). Distance to water was then calculated for each pixel as the Euclidean distance to the nearest body of water using the distance command in the R package raster (48). We extracted the mean value of each of the above variables for each protected area using zonal statistics.

Distance to protected areas and protected area size.

We rasterized the protected areas of Costa Rica (https://www.protectedplanet.net/) with those pixels being composed of at least 50% protected area at a 915 m resolution being assigned a value of 1. We then calculated the distance to the nearest protected area using the “distance” command in the R package raster (48). The total surface area for the protected areas, which we downloaded from the world’s protected area network database, was also used as a predictor variable ().

Biodiversity variables.

We downloaded occurrence points for terrestrial vertebrate species (mammals, reptiles, amphibians, and birds) in Costa Rica from the Global Biodiversity Information Facility (GBIF) in July 2021 (GBIF.org, https://doi.org/10.15468/dl.4a3zpk). We removed records that had geospatial issues as noted by GBIF (outside of assumed datum or rounded coordinates). We eliminated all occurrence records not resolved to species, with spatial resolution less than 1 km. We also excluded any species for which we had fewer than 25 occurrence records. We only included observations from 2005 to 2015 to match the temporal extent of the Flickr and eBird data (see below).

With these data, we estimated spatial terrestrial vertebrate species richness in Costa Rica using MaxEnt spatial modeling software version 3.1 (50), which generates probabilistic range maps using the climatic and land cover conditions present where species are found. Models were fit using four climatic variables from WorldClim (annual temperature bio1, temperature seasonality bio4, annual precipitation bio12, and precipitation seasonality bio15). Climate data were acquired at a 30 arc-second resolution across the entire country (46). These files were reprojected to UTM zone 16N from their native projection system (Wideband Global SATCOM, 1984). All spatial products were also reprojected to the UTM zone 16N system using bilinear interpretation.

We gathered Landsat satellite images (United States Geological Survey, United States) to create a land-use/land-cover map. We used fully-constrained spectral mixture analysis to create subpixel fractional cover maps of soil (fcover1), photosynthetic vegetation (fcover2), and impervious surfaces (fcover3), also from Landsat (75). Soil and impervious spectra were resampled from a global spectral library (76), and vegetation spectra were simulated using PROSAIL (77). We created percent tree cover maps using random forest regression, with high resolution tree cover maps from (78) as training data and Landsat/fcover as covariates. Each layer was created at 30 m resolution and upscaled them to 915 m resolution (47). Thus, models were also fit using the land cover maps. Open source code to create these maps is available in GitHub (https://github.com/earth-chris/earthlib/).

To account for uneven spatial distribution of sampling effort, we created a bias raster using all GBIF occurrences within the extent of Costa Rica. After applying the same quality controls described above, we reprojected the coordinates to UTM zone 16N. We then used the "kde2d" function of the MASS package with a Gaussian decay function and a 5-km SD to create a raster of the relative density of GBIF points across the landscape. We used this map to randomly sample 10,000 pseudo absences from the background, with the likelihood of a pixel being chosen being weighted by the value of that pixel in the bias map.

We allowed the model to be fit with linear, quadratic, multiplicative, and hinge terms using standard penalization techniques for more complex terms. The models were fit with four bioclimatic variables (mean annual temperature, mean annual precipitation, temperature seasonality, and precipitation seasonality), percent soil cover, percent photosynthetic cover, and percent tree cover. The run was replicated five times for each species, withholding 20% of the data for testing the models fit using AUC. All other parameters were kept as the MaxEnt default options. We excluded any species where the AUC of the species distribution model when evaluating the test dataset was less than 0.75, indicating a poor model fit. A total of 765 vertebrate species were ultimately included in the analysis (birds = 699 spp, mammals= 23 spp, reptiles = 21 spp, amphibians= 22 spp, see ). All models were conducted using the dismo package in R (51). To subset species that were classified as threatened and endemic birds, we took the list of birds (n = 699) and consulted the BirdLife International and International Union for Conservation of Nature (IUCN) red list status. We classified a species as threatened or endemic if the species or a subspecies was endemic to Costa Rica, endemic to Central America or listed as vulnerable and endangered by IUCN. In total, we found 279 species that belonged in this category (). We used these three biodiversity variables (1-bird species richness, 2-richness of threatened and endemic birds, and 3-amphibian, mammals, reptiles species richness) in separate models to avoid collinearity given the high correlation between them ().

We extracted the mean prevalence for each species, by first making a binary presence/absence map of each species using a thresholding method common in the literature which sets the probability of a false-positive and a false-negative to be equal (79). These were then summed to create species richness maps using species belonging to each biodiversity variable (either bird species richness, endemic and threatened birds, or richness of amphibians, reptiles, and mammals).

Tourism.

We measured tourism via Flickr photographs. Flickr is a social media platform for sharing photographs, some of which are tagged with the date and location where the user captured the image. Visitation was estimated by summing unique Flickr users who shared at least one image from a location per day. This metric, known as photograph-user-days (PUD) (19) corrects for potential biases created by Flickr users who upload multiple photographs from one location per day (19), and it is the most common measure of visitation for research on recreation ecosystem services (79). We used the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) visitation model to compute PUD across all the available years of data (2005–2017) (80). We used two separate InVEST runs to generate PUDs for each protected area and for each 915 m × 915 m pixel country-wide. We also created a threshold of 4 km2 as the minimum protected area size. We did so because we identified that protected areas smaller than 4 km2 had zero or very few Flickr photograph uploads and were inducing skewness in the data via zero-inflated patterns.

To measure birdwatching tourism, we used all eBird observations made from 2005 to 2017 (81). eBird is an online platform where birdwatchers share geolocated checklists of birds observed at a specific time. As we did for Flickr PUDs, we only included unique eBird checklist-user-days (CUD), so that a given user contributed no more than a single count per day in an individual area. We summed CUD across all years in each protected area for the protected area analysis and per pixel (915 m × 915 m) for the country-wide analysis.

Data Analysis.

We used MaxEnt software (82) to map how tourism can be explained by predictor variables. We used species distribution models and instead of species occurrence variables we used the two tourism responses: Flickr PUD and eBird CUD. We applied square root transformations to distance variables (distance to water and distance to protected areas), and applied a log transformation to hotel density. We did so to conform with model assumptions and to assess goodness of fit for linear regressions (83) (). We ran ten models in total, out of which five were at the country-wide scale and five at the protected area scale (see ).

We fit models using distance to roads, hotel density, distance to water, protected area size (for models within protected areas) or distance to protected areas (for the country-wide analyses), and the three biodiversity variables independently (see for all ten models). For all models, the selected five layers used as predictor variables did not show collinearity (variance inflation factor [VIF] <2, ) as measured by the R package usdm. For the models that considered tourism within protected areas, we clipped all input layers using polygons provided by the Protected Areas of Costa Rica. Unlike our species distribution models for biodiversity, we constrained the model to only select linear and/or quadratic terms. Models were replicated ten times, withholding 10% of the data for validation in a k-fold pattern. Pseudo absences were chosen at random across the entirety of the extent (either the whole country for country-wide models or just within protected areas). Otherwise default options were used in the models. All models were conducted using the R dismo package (84).

In the main text, we report results of all ten models. In the , we report summary statistics for each model (e.g., model AUC, percent contribution, permutation importance for each variable, see ) and show graphically how each variable affects tourism (Figs. 1 and 2). We show this by plotting on the x-axis each predictor variable and on the y-axis the expected tourism value for MaxEnt models constructed using only this variable.

We projected tourism as either Flickr PUD or eBird CUD across the entirety of the country using the predict function in the dismo package in R (Fig. 3 ). Lastly, we created maps of the difference between the models created with all variables and those which withheld either biodiversity (bird richness, threatened and endemic birds, or amphibians, reptiles, and mammals) (Fig. 3 ), or infrastructure variables (distance to roads and hotel density) (Fig. 3 ) to see how excluding these variables affected the predicted patterns of tourism.