Potential Current and Future Habitats of Capra ibex in the European Alps: Does Capra ibex really profit from Climate Change?

Potential Current and Future Habitats of Capra ibex in the European Alps: Does Capra ibex really profit from Climate Change?


Samuel Hoffmann

Graduate Program (M.Sc.) Global Change Ecology

within

Module B5 – Global Change Impacts on Species Distribution

31.07.2014


These analysis have been conducted within a semester course to learn how to apply different techniques in ecology. The codes can be re-used, but please keep in mind that all analysis have been done within a student project – no guarantee of completeness or correctness can be given.


Introduction

Climate induced environmental change leads to loss of biodiversity, particularly in rather extreme ecosystems such as the alpine ecosystems in the European Alps (Pauli et al., 2012). Major reasons are the disruption of seasonal synchronization of interdependent ecological processes and trophic interactions, even though abiotic drivers might have been the trigger, but distinguishing biotic from abiotic causes is difficult, since highly spatial and temporal resolute species data is sparse (Post et al., 2008, Büntgen et al., 2014).

In contrast, seasonal plant growth will be enhanced due to rising springtime temperatures, and earlier and faster snow melt predicted by climate models. This would, in turn, improve forage value of habitat vegetation for many mammal species including Capra ibex (fig. 1), which is endemic to Europe and native to the European Alps. Hence, it is assumed that alpine ibex is among the winners of climate change with regards to a study of Büntgen et al. (2014), comparing annual climate data and horn growth as an indicator for individual performance. Nevertheless, quantity and quality of food resources are just few of many factors influencing population viability, and are, moreover, not always correlated with species growth and survival (Pettorelli et al., 2007). Even if the species is not considered to be threatened nowadays, there is a certain risk of extinction due to low genetic variability, the founder effect and minimum viable populations (Shackelton, 1997, Maudet et al., 2002). Inbreeding depression, hybridization, habitat fragmentation, epizootics, human disturbances such as tourism and poaching are major threats to Capra ibex (Aulagnier et al., 2008).

Fig. 1: Capra ibex in its characteristic environment. Source: Andres Eraso-Keller (https://www.flickr.com/photos/59558842@N05/).

Habitat availability and suitability according to size, spatial structure, and climate, are other important aspects for species vitality in terms of population density, resource availability, reproduction rate and breeding ground (Amon, 1956, Wiersema, 1989, Grignolio et al., 2003, Grignolio et al., 2004, Jacobson et al., 2004, Aulagnier et al., 2008). Therefore, this investigation concerns modelling of potential actual and future habitats for Capra ibex in the European Alps, implementing topographic and climatic data, in order to depict variation in habitat suitability and availability owing to climate change. Does Capra ibex really profit from climate change in the long run as suggested by Büntgen et al. (2014)? Because, intuitively I would assume that climate change leads to a severe reduction of suitable habitats, since at some point an upper altitudinal boundary in alpine ecosystems hinders alpine ibex migrating further up, in order to compensate changing temperature and precipitation patterns with altitude. Bringing this together with already known threats as mentioned above, Capra ibex could be prone to extinction in future again.

Methods

In the European Alps Capra ibex was prone to extinction in the 19th century. About 100 individuals survived in “Gran Paradiso”, former private land for hunting purposes owned by King Vittorio Emanuele III. and these days a National Park in north-western Italy (Willmann, 2014). A breeding program, starting at the beginning of the 20th century, facilitated re-introduction and re-settlement of Capra ibex in almost the entire Alps arc. The IUCN Red List of Threatened Species contains information about Capra ibex, which currently describes 151 areas of occurrence and c. 31670 individuals in the Alps. Out of the 151 species occurrence polygons, 31670 species occurrence points were randomly selected for further analyses. At first, literature research provided information about meaningful environmental variables for alpine ibex such as snow depth or respectively winter precipitation and temperature, temperature and precipitation in spring and early summer, insolation (i.e. topographic slope and aspect) and land cover (Amon, 1956, Wiersema, 1989, Grignolio et al., 2003, Grignolio et al., 2004, Jacobson et al., 2004, Aulagnier et al., 2008, Grotan et al., 2008,). In addition, distances to broad roads and large urban areas were incorporated, which are commonly considered as disturbances for many mammal species (EEA & FOEN, 2011). Afterwards freely available environmental data, which corresponds to these variables, was searched. Climate data was derived from “WorldClim – Global Climate Data” (Hijmans et al., 2005). Future climate data was predicted by the HadGEM2_AO General Circulation Model for the year 2070, taking the Representative Concentration Pathway 4.5. Thereby six so-called “bioclimatic” variables were selected: Mean temperature of wettest quarter (BIO 8), mean temperature of driest quarter (BIO 9), mean temperature of warmest quarter (BIO 10), mean temperature of coldest quarter (BIO 11), precipitation of warmest quarter (BIO 18), precipitation of coldest quarter (BIO 19). Topographical data was derived from a digital elevation model, given by the International Centre for Tropical Agriculture (CIAT) (Jarvis et al., 2008), and land cover data was retrieved from the “Global Land Cover 2000” database (European Commission Joint Research Centre, 2003). Data for urban areas, rivers and lakes stems from “Mapcruzin.com” (Meuser, 2014). Road data was obtained from “DIVA-GIS” database (Hijmans, 2014), boundary data for the European Alps from the European Soil Portal (European Commission, Joint Research Centre, 2012). Spatial data resolution is or was fitted to 1 km. Utilized data and related webpages are listed in table 1.

Tab. 1: Utilized data and related webpages.
Data Webpage
Species Occurrence http://www.iucnredlist.org
Climate http://www.worldclim.org
Digital Elevation Model http://www.srtm.csi.cgiar.org
Land cover http://www.bioval.jrc.ec.europa.eu
Urban Areas, Rivers, Lakes http://www.mapcruizin.com
Roads http://www.diva-gis.org/
Boarder http://www.eusoils.jrc.ec.europa.eu

Collinearity was tested and removed. Variables were excluded according to a correlation value of >0.7 (see Appendix). The parameters of mean temperature of warmest quarter (BIO 10), precipitation of coldest quarter (BIO 19), slope, aspect, landcover and distances to roads and urban areas are remaining. Three different modelling algorithms were applied and evaluated: General Additive Model, Maxent Model, and Random Forest Model. All analyses were conducted using open source software “R Studio” and “GRASS GIS”. In order to comprehend and reproduce programming procedure, the used R script is attached. For further details see Appendix.

Results

A map, showing the probability of species occurrence, calculated for each of the three models as well as for current and future climate conditions, was produced (fig. 2). Furthermore, the probability differences of each model between presence and future were mapped. Supplementary, the mean and standard deviation of all three model outputs for presence, future and the differences between are illustrated. That means, on the one hand side, that more reliable results are generated, since they are based on the output of all three models. On the other hand side, the map of standard deviations is supposed to identify areas of inconsistent model outcomes. However, one has to be aware of the fact that the computed probability values of the different models cannot be weighted in the same way, because the model algorithms act differently. Nonetheless, mapping mean and standard deviation gives an impression, where areas of more significant results can be found, even though absolute probability values should not be taken for granted.

plot of chunk unnamed-chunk-2

Fig. 2: Maps showing the probability of species occurence, calculated by a General Additive Model (GAM, C(ROC)= 0.95), Random Forest (C(ROC)= 0.95) and Maxent Model (C(ROC)= 0.94) for actual and future climate (2070, HadGEM2-AO GCM, RCP 4.5), as well as the differences between both climatic conditions (Future-Actual). Additionally, the mean and standard deviation (Sd) of all three models for presence, future and differences between are illustrated. Actual species occurences are marked by black dots. The black line reflects the boarder of the European Alps.

It can be seen that there is a general trend towards a reduction of species occurrence probability from presence to future climate conditions, which conforms to a reduction of overall habitat size. None of the models predict an increase in occurrence probability. Especially low valleys are affected. Still describing rather suitable habitats at present, future climate scenario might convert valleys to rather unlikely habitats. The GAM model indicates greatest loss of areas with high probability of species occurrence, followed by the Random Forest model and the Maxent model. According to the C-indices or areas under the receiver operating characteristic curve model results are reliable (tab.2).

Tab. 2: Applied Models and their C-indices or areas under the receiver operating characteristic curve, respectively.
Model C(ROC)
GAM 0.95
Random Forest 0.95
Maxent 0.94

Highest probabilities of c. 100 % are concentrated at the Alpine divide, the highest part of the Alps, which are also denoted by highest values of standard deviation. In fact, standard deviations remain relatively low, while means remain relatively high, so that one can rely on the model output. Yet, a certain model discordance cannot be neglected, mainly caused by the models' different approaches. Interestingly the standard deviations of the model differences are highest in areas mostly in front of the Alpine divide and mountainous regions outside the alps. This variation in the model results might also stem from the diverse model algorithms, which choose different environmental input factors to be most representative for the species existence.

Hence, it is necessary to assess the relative importance of each environmental parameter, influencing the model output (see Appendix). As the occurrence probability patterns go along with altitudinal patterns, and altitude and temperature are strongly correlated, it seems that temperature is the main environmental driver in the models, which is confirmed by model statistics (see Appendix). “Mean temperature of warmest quarter” (BIO 10) is the most important factor in the Random Forest and Maxent model. In summary, apart from temperature related parameters, slope also matters a lot, whereas the bio-climatic variable concerning precipitation in coldest quarter (BIO 19), distances to roads and urban areas as well as aspect are of much less importance. Only in the Maxent model land cover is supposed to be relevant.

Another obvious aspect is that the species probably won't be able to survive in most of the lower mountainous regions outside the European Alps in 2070. A potential introduction of the species to these regions would be pointless for the purpose of nature conservation. Moreover, some recent areas, where Capra ibex appears, are expected to become inadequate for species' continuance. For this reason, revised and appropriate wildlife management strategies should be developed to ensure the species survival in these areas.

Discussion

The results suggest that the current area of suitable habitats of Capra ibex in the European Alps is clearly decreasing under future climate conditions. And yet the selected climate scenario is predicting moderate rates of change compared to other scenarios. In fact, the abandonment of traditional agriculture in the Alps leads to even more habitat loss, because of vegetation succession from pasture to forest. For this reason, even if rising spring time temperatures enhance population viability, the existence of Capra ibex might be threatened by changing habitat conditions at some point in space and time. Apart from that, it is still a subject under discussion how alpine vegetation and, accordingly, food resources will change and affect alpine ungulates. In any case, the species will have to react to survive, either in terms of flexibility, i.e. phenotypic change or habitat change, acclimatization, learning processes, or with respect to genetic adaptation, which probably takes too long contrasted to the rate of climate change and the given bottleneck effect within the species' gene pool. Specific management strategies could support the species to accomplish that. Less hunting activities, feeding during winter, continuing traditional agriculture, decreasing human induced disturbances, reducing habitat fragmentation by land use changes, preventing hybridization by excluding domestic sheep and goat in re-introduction areas, or proceeding reintroduction by considering suitable, future habitats and genetic diversity are just some starting points one can think about (Aulagnier et al, 2008). Presenting sophisticated management plans, however, would exceed the purpose of this study.

Although many investigations state that snow depth is an important environmental parameter for ibex' vitality (Wiersema, 1989, Grignolio et al., 2003, Grignolio et al., 2004, Jacobson et al., 2004, Aulagnier et al., 2008, Grotan et al., 2008,), the models of this study do not. It points out that temperature and slope are generally associated with species occurrence at most. A reason for that disagreement might be that sampled absence points are missing. Especially the mean temperature in the warmest quarter of the year is highly rated, which is surprising, because summer temperatures are hardly mentioned to be meaningful in literature concerning alpine ibex. Maybe it is not even temperature directly, which drives species occurrence, but the corresponding altitude, which influences the type of ecosystem, which appears there and on which alpine ibex depends, in turn. Further errors and inaccuracy is certainly inherent in the method of randomly choosing occurrence points out of the species polygon data, since spatial data resolution of 1 km is simply too coarse to reflect within-habitat variability in reality. Thus, land cover classes might matter less than assumed before. Besides, Capra ibex changes habitats throughout the year, which is not accounted for in these models. In winter it prefers lower elevations including open rocks, pastures and forest patches, whereas in summer it moves above the tree line (Aulagnier et al., 2008). Also collinearity, which was not completely removed, has an effect on the result, which was not examined.

Continuative analyses should incorporate that. Other, maybe even less variables for more significant outcomes and better evaluation can be tested as well. Resolution should be increased and other model algorithms can be applied. The consequence of a different method to select species points can be reviewed, too.

Conclusion

Finally, claiming that Capra ibex is a winner under climate change in the European Alps is surely exaggerated, since horn growth of some individuals from few populations might be a proxy for individual performance under short-time climate variations, but not a appropriate proxy for estimating the species' entire potential to adapt to diverse future environmental changes in the long run, which are triggered by climate change. In this regard, this study has shown that habitats of the alpine ibex are considerably decreasing in future, even though the methodology is not yet fully developed. Hence, the results affirm that Capra ibex will undergo diverse positive and negative environmental changes, whose total impact on the species is still difficult to assess. Further research is desirable. Because, if we stop paying attention to such impressive species, we might watch another species becoming extinct in the wild quietly.

Acknowledgements

At first, I want to thank Martin Wegmann and Björn Reineking for their expertise and support in developing this study. Also many thanks to my colleague Joseph Premier for inspiring and stimulating discussions. This investigation was conducted within course “B5 – Global Change Impacts on Species Distribution” in the Graduate Program (M. Sc.) “Global Change Ecology” within the Elite Network of Bavaria at the University of Bayreuth.

References

  • Amon R. (1956): Klimatische Grenzen einer Wiederansiedlung des Alpensteinbockes,Capra i. ibex, Linné 1758, in den Ostalpen, Zeitschrift für Jagdwissenschaft 5(3):132-137.

  • Aulagnier S., Kranz A., Lovari S., Jdeidi T., Masseti M., Nader I., de Smet, K., Cuzin F. (2008): Capra ibex, In: IUCN 2013, IUCN Red List of Threatened Species, Version 2013.2, http://iucnredlist.org, Downloaded on 10 June 2014.

  • Büntgen U., Liebhold A., Jenny H., Mysterud A., Egli1 S., Nievergelt D., Stenseth N. C., Bollmann K. (2014): European springtime temperature synchronises ibex horn growth across the eastern Swiss Alps, Ecology Letters 17(3):303-313.

  • European Environment Agency (EEA) & Swiss Federal Office for the Environment (FOEN) (2011): Landscape fragmentation in Europe, Joint EEA-FOEN report, Schultz Grafisk, Copenhagen.

  • European Commission, Joint Research Centre (2012): European Soil Portal database, http://www.eusoils.jrc.ec.europa.eu

  • European Commission, Joint Research Centre (2003): Global Land Cover 2000 database, http://bioval.jrc.ec.europa.eu/products/glc2000/glc2000.php

  • Grignolio S., Parrini F., Bassano B., Luccarini S., Apollonio M. (2003): Habitat selection in adult males of Alpine ibex, Capra ibex ibex, Folia Zoologica 52(2):113-120.

  • Grignolio S., Rossi I., Bassano B., Parrini F., Apollonio M. (2004): Seasonal variations of spatial behaviour in female Alpine ibex (Capra ibex ibex) in relation to climatic conditions and age, Ethology Ecology and Evolution 16:255-264.

  • Grotan V., Saether B.-E., Filli F., Engen S. (2008): Effects of climate on population fluctuations of ibex, Global Change Biology 14(2):218-228.

  • Hijmans, R.J., S.E. Cameron, J.L. Parra, P.G. Jones and A. Jarvis (2005): Very high resolution interpolated climate surfaces for global land areas, International Journal of Climatology 25:1965-1978.

  • Hijmans, R.J. (2014): http://www.diva-gis.org/

  • Hirzel A. H., Hausser J., Chessel D., Perrin N. (2002): Ecological-Niche Factor Analysis: How to compute Habitat-Suitability Maps without Absence Data? Ecology 83(7):2027-2036.

  • Jacobson A. R., Provenzale A., von Hardenberg A., Bassano B., Festa-Bianchet M. (2004): Climate Forcing and Density Dependence in a Mountain Ungulate Population, Ecology 85(6):1598-1610.

  • Jarvis A., H.I. Reuter, A. Nelson, E. Guevara (2008): Hole-filled seamless SRTM
    data V4, International Centre for Tropical Agriculture (CIAT), available from
    http://srtm.csi.cgiar.org

  • Maudet, C., Miller, C., Bassano, B., Breitenmoser-Würsten, C., Gauthier, D., Obexer-Ruff, G., Michallet, J., Taberlet, P. and Luikart, G. (2002): Microsatellite DNA and recent statistical methods in wildlife conservation management: applications in Alpine ibex [Capra ibex (ibex)], Molecular Ecology 11(3):421-436.

  • Meuser, M. (2014): http://www.mapcruzin.com

  • Mignatti A., Casagrandi R., Provenzale A., von Hardenberg A., Gatto M. (2012): Sex- and age-structured models for Alpine ibex Capra ibex ibex population dynamics, Wildlife Biology 18(3):318-332.

  • Pauli, H., Gottfried, M., Dullinger, S., Abdaladze, O., Akhalkatsi, M., Alonso, J.L.B. et al. (2012): Recent plant diversity changes on Europe's mountain summits, Science 336:353-355.

  • Pettorelli N., Pelletier F., von Hardenberg A., Festa-Bianchet M., Côté S. D. (2007): Early Onset of Vegetation Growth vs. Rapid Green-Up: Impacts on Juvenile Mountain Ungulates, Ecology 88(2):381-390.

  • Post, E., Pedersen, C., Wilmers, C.C. & Forchhammer, M.C. (2008): Warming, plant phenology and the spatial dimension of trophic mismatch for large herbivores, Proceedings of the Royal Society B 27:2005-2013.

  • Shackleton, D.M. (1997): Wild Sheep and Goats and Their Relatives: Status Survey and Conservation Action Plan for Caprinae. In: D. M. Shackleton (ed.), Wild sheep and goats and their relatives. Status survey and conservation action plan for Caprinae, IUCN/SSC Caprinae Specialist Group, Gland, Switzerland and Cambridge, UK.

  • Wiersema G. (1989): Climate and Vegetation Characteristics of Ibex Habitats in the European Alps, Mountain Research and Development 9(2):119-128.

  • Willmann, U. (2014): Steinbock: Protz im Gebirge, Die Zeit, Nr. 25, http://www.zeit.de/2014/25/steinbock-alpen-klimaerwaermung

Appendix

Set working directory.

Load required packages:

library(rms) 
library(raster)
library(mgcv)
library(gam)
library(randomForest)
library(dismo)
library(rgdal)
library(ellipse)
library(rJava)
library(XML)
library(sp)
library(proj4)
library(rgeos)    
library(maptools)
library(GISTools)
library(rasterVis)

Load function “varImpBiomod.R” within working directory:

source("varImpBiomod.R") 

Load study area boarder and change projection:

alps_boundary <- readOGR("./GIS/urban_areas_alps","AlpineConvention") 
alps_boundary <- spTransform(alps_boundary, CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

Load environmental variables, which are ecologically meaningful for Capra ibex (slope, aspect, land cover):

alps_dem1 <- raster("./GIS/CSI_dem/srtm_36_03.tif")  
alps_dem2 <- raster("./GIS/CSI_dem/srtm_36_04.tif")   
alps_dem3 <- raster("./GIS/CSI_dem/srtm_37_03.tif")   
alps_dem4 <- raster("./GIS/CSI_dem/srtm_37_04.tif")   
alps_dem5 <- raster("./GIS/CSI_dem/srtm_38_03.tif")   
alps_dem6 <- raster("./GIS/CSI_dem/srtm_38_04.tif")   
alps_dem7 <- raster("./GIS/CSI_dem/srtm_39_03.tif")   
alps_dem8 <- raster("./GIS/CSI_dem/srtm_39_04.tif")   
alps_dem9 <- raster("./GIS/CSI_dem/srtm_40_03.tif")   
alps_dem10 <- raster("./GIS/CSI_dem/srtm_40_04.tif")   

alps_dem <- merge(alps_dem1, alps_dem2, alps_dem3,alps_dem4, alps_dem5, alps_dem6,alps_dem7, alps_dem8, alps_dem9, alps_dem10)

plot(alps_dem)
plot(alps_boundary, add=T)

alps_dem90 <- crop(alps_dem,extent(alps_boundary)+1)

plot(alps_dem90)
plot(alps_boundary, add=T)

writeRaster(alps_dem90,filename="./GIS/alps_dem90.grd")
alps_dem90 <- stack("./GIS/alps_dem90.grd")

alps_slope <- terrain(alps_dem90,opt="slope",unit="degrees",neighbors=8)
alps_aspect <- terrain(alps_dem90,opt="aspect",unit="degrees")

plot(alps_slope)
plot(alps_aspect)

extent(alps_slope)
extent(alps_aspect)

lacov <- raster("./GIS/GRID/glc_eu_v2/hdr.adf")
projection(lacov)
plot(lacov)
alps_lacov <- crop(lacov,extent(alps_boundary)+1)
plot(alps_lacov)
alps_dem90 <- raster("./GIS/alps_dem90.grd")

Upscale slope and aspect from 90 m to 1 km (resolution of climate data):

slope_aggregate <- aggregate(alps_slope, fact=7, fun=mean)
extent(slope_aggregate)
slope_resample <- resample(slope_aggregate, alps_bioclim)

aspect_aggregate <- aggregate(alps_aspect, fact=7, fun=mean)
extent(aspect_aggregate)
aspect_resample <- resample(aspect_aggregate, alps_bioclim)

Load data for cities, lakes, rivers in the Alps including a buffer of approx. 1 km for urban areas and roads:

alps_urban <- readOGR("./GIS/urban_areas_alps","urban_alps") 
alps_urban <- spTransform(alps_urban, CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
plot(alps_urban)
class(alps_urban)
alps_urban <- rasterize(alps_urban, slope_aggregate)
alps_urban[is.na(alps_urban)]<-0
alps_urban
writeRaster(alps_urban,filename="./GIS/alps_urban.grd", overwrite=T)

alps_urban_buf <- gBuffer(alps_urban, width=0.01)
plot(alps_urban_buf)
alps_urban_buf <- rasterize(alps_urban_buf, slope_aggregate)
alps_urban_buf[is.na(alps_urban_buf)]<-0
alps_urban_buf
writeRaster(alps_urban_buf,filename="./GIS/alps_urban_buf.grd", overwrite=T)

alps_lakes <- readOGR("./GIS/lakes_alps","lake_alps") 
alps_lakes <- spTransform(alps_lakes, CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
plot(alps_lakes)
class(alps_lakes)
alps_lakes <- rasterize(alps_lakes, slope_aggregate)
alps_lakes[is.na(alps_lakes)]<-0
alps_lakes
writeRaster(alps_lakes,filename="./GIS/alps_lakes.grd", overwrite=T)

alps_rivers <- readOGR("./GIS/rivers_alps","river_alps") 
alps_rivers <- spTransform(alps_rivers, CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
plot(alps_rivers)
class(alps_rivers)
alps_rivers <- rasterize(alps_rivers, slope_aggregate)
alps_rivers[is.na(alps_rivers)]<-0
alps_rivers
writeRaster(alps_rivers,filename="./GIS/alps_rivers.grd", overwrite=T)

#alps_roads <- readOGR("./GIS/alps-roads-shape","roads") 
#alps_roads <- spTransform(alps_roads, CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
#plot(alps_roads)
#class(alps_roads)
#alps_roads <- rasterize(alps_roads, slope_aggregate)
#alps_roads[is.na(alps_roads)]<-0
#alps_roads
#writeRaster(alps_roads,filename="./GIS/alps_roads.grd", overwrite=T)

alps_roads1 <- readOGR("./GIS/DEU_rds","DEU_roads") 
alps_roads2 <- readOGR("./GIS/AUT_rds","AUT_roads") 
alps_roads3 <- readOGR("./GIS/CHE_rds","CHE_roads") 
alps_roads4 <- readOGR("./GIS/FRA_rds","FRA_roads") 
alps_roads5 <- readOGR("./GIS/HUN_rds","HUN_roads") 
alps_roads6 <- readOGR("./GIS/ITA_rds","ITA_roads") 
alps_roads7 <- readOGR("./GIS/LIE_rds","LIE_roads") 
alps_roads8 <- readOGR("./GIS/MCO_rds","MCO_roads") 
alps_roads9 <- readOGR("./GIS/SVN_rds","SVN_roads") 
alps_roads10 <- readOGR("./GIS/SVK_rds","SVK_roads") 
alps_roads11 <- readOGR("./GIS/CZE_rds","CZE_roads") 
alps_roads12 <- readOGR("./GIS/HRV_rds","HRV_roads") 
alps_roads13 <- readOGR("./GIS/BIH_rds","BIH_roads") 

alps_roads1 <- gBuffer(alps_roads1, width=0.01)
alps_roads2 <- gBuffer(alps_roads2, width=0.01)
alps_roads3 <- gBuffer(alps_roads3, width=0.01)
alps_roads4 <- gBuffer(alps_roads4, width=0.01)
alps_roads5 <- gBuffer(alps_roads5, width=0.01)
alps_roads6 <- gBuffer(alps_roads6, width=0.01)
alps_roads7 <- gBuffer(alps_roads7, width=0.01)
alps_roads8 <- gBuffer(alps_roads8, width=0.01)
alps_roads9 <- gBuffer(alps_roads9, width=0.01)
alps_roads10 <- gBuffer(alps_roads10, width=0.01)
alps_roads11 <- gBuffer(alps_roads11, width=0.01)
alps_roads12 <- gBuffer(alps_roads12, width=0.01)
alps_roads13 <- gBuffer(alps_roads13, width=0.01)

alps_roads1 <- rasterize(alps_roads1, slope_aggregate)
alps_roads2 <- rasterize(alps_roads2, slope_aggregate)
alps_roads3 <- rasterize(alps_roads3, slope_aggregate)
alps_roads4 <- rasterize(alps_roads4, slope_aggregate)
alps_roads5 <- rasterize(alps_roads5, slope_aggregate)
alps_roads6 <- rasterize(alps_roads6, slope_aggregate)
alps_roads7 <- rasterize(alps_roads7, slope_aggregate)
alps_roads8 <- rasterize(alps_roads8, slope_aggregate)
alps_roads9 <- rasterize(alps_roads9, slope_aggregate)
alps_roads10 <- rasterize(alps_roads10, slope_aggregate)
alps_roads11 <- rasterize(alps_roads11, slope_aggregate)
alps_roads12 <- rasterize(alps_roads12, slope_aggregate)
alps_roads13 <- rasterize(alps_roads13, slope_aggregate)

alps_roads_buf <- mosaic(alps_roads1, alps_roads2, alps_roads3, alps_roads4, alps_roads5, alps_roads6, alps_roads7, alps_roads8, alps_roads9,alps_roads10,alps_roads11,alps_roads12,alps_roads13, fun=mean)
alps_roads_buf <- crop(alps_roads_buf,extent(alps_boundary)+1)
alps_roads_buf <- resample(alps_roads_buf, slope_aggregate)
alps_roads_buf[is.na(alps_roads_buf)]<-0
alps_roads_buf

alps_roads <- mosaic(alps_roads1, alps_roads2, alps_roads3, alps_roads4, alps_roads5, alps_roads6, alps_roads7, alps_roads8, alps_roads9,alps_roads10,alps_roads11,alps_roads12,alps_roads13, fun=mean)
alps_roads <- crop(alps_roads,extent(alps_boundary)+1)
alps_roads <- resample(alps_roads, slope_aggregate)
alps_roads[is.na(alps_roads)]<-0
alps_roads

plot(alps_roads_buf)
plot(alps_boundary, add=T)

writeRaster(alps_roads_buf,filename="./GIS/alps_roads_buf.grd", overwrite=T)
alps_roads_buf <- stack("./GIS/alps_roads_buf.grd")

writeRaster(alps_roads,filename="./GIS/alps_roads.grd", overwrite=T)
alps_roads <- stack("./GIS/alps_roads.grd")

alps_ulrr_buf <- overlay(alps_urban_buf, alps_lakes, alps_rivers, alps_roads_buf, fun=sum)
plot(alps_ulrr_buf)
alps_ulrr_buf[which(alps_ulrr_buf[]>0)]<-1
writeRaster(alps_ulrr_buf,filename="./GIS/alps_ulrr_buf.grd", overwrite=T)

alps_ulrr <- overlay(alps_urban, alps_lakes, alps_rivers, alps_roads, fun=sum)
plot(alps_ulrr)
alps_ulrr[which(alps_ulrr[]>0)]<-1
writeRaster(alps_ulrr,filename="./GIS/alps_ulrr.grd", overwrite=T)

Calculate distances from environmental data to roads and urban areas:

alps_ur <- overlay(alps_urban, alps_roads, fun=sum)
plot(alps_ur)
alps_ur[which(alps_ur[]>0)]<-1
writeRaster(alps_ur,filename="./GIS/alps_ur.grd", overwrite=T)
alps_ur <- raster("./GIS/alps_ur.grd")
writeRaster(alps_ur,filename="./GIS/alps_ur.tif", overwrite=T)
alps_ur <- raster("./GIS/alps_ur.tif")

alps_ur_dist <- alps_ur
alps_ur_dist[which(alps_ur_dist[]==0)]<-NA
alps_ur_dist[alps_ur_dist[]==0]<-NA

plot(alps_ur_dist) 
str(alps_ur_dist)
is.na(alps_ur_dist)
?is.na

dist <- raster::distance(alps_ur_dist)

dist <- raster("./GIS/distance.tif")
projection(dist) <- projection(alps_ulrr)
class(dist)
plot(dist)
writeRaster(dist, filename="./GIS/alps_dist.tif", overwrite=T)
alps_dist <- raster("./GIS/alps_dist.tif")

Load data again:

alps_ulrr <- stack("./GIS/alps_ulrr.grd")
alps_ulrr_buf <- stack("./GIS/alps_ulrr_buf.grd")

Load climate data, which is ecologically meaningful for Capra ibex (according to literature: snow depth, summer temperature, spring temperature, precipitation in early summer, precipitation from Dec-Jun, insolation etc.):

bioclim <- raster::getData("worldclim", var = "bio", res = 0.5, lon=10, lat=45)
bioclim <- subset(bioclim, c("bio8_16", "bio9_16", "bio10_16", "bio11_16","bio18_16","bio19_16"))

plot(bioclim,1)
plot(alps_boundary,add=T)

bioclim2070_1 <- raster("./GIS/bioclim2070/hd45bi701.tif")
bioclim2070_2 <- raster("./GIS/bioclim2070/hd45bi702.tif")
bioclim2070_3 <- raster("./GIS/bioclim2070/hd45bi703.tif")
bioclim2070_4 <- raster("./GIS/bioclim2070/hd45bi704.tif")
bioclim2070_5 <- raster("./GIS/bioclim2070/hd45bi705.tif")
bioclim2070_6 <- raster("./GIS/bioclim2070/hd45bi706.tif")
bioclim2070_7 <- raster("./GIS/bioclim2070/hd45bi707.tif")
bioclim2070_8 <- raster("./GIS/bioclim2070/hd45bi708.tif")
bioclim2070_9 <- raster("./GIS/bioclim2070/hd45bi709.tif")
bioclim2070_10 <- raster("./GIS/bioclim2070/hd45bi7010.tif")
bioclim2070_11 <- raster("./GIS/bioclim2070/hd45bi7011.tif")
bioclim2070_12 <- raster("./GIS/bioclim2070/hd45bi7012.tif")
bioclim2070_13 <- raster("./GIS/bioclim2070/hd45bi7013.tif")
bioclim2070_14 <- raster("./GIS/bioclim2070/hd45bi7014.tif")
bioclim2070_15 <- raster("./GIS/bioclim2070/hd45bi7015.tif")
bioclim2070_16 <- raster("./GIS/bioclim2070/hd45bi7016.tif")
bioclim2070_17 <- raster("./GIS/bioclim2070/hd45bi7017.tif")
bioclim2070_18 <- raster("./GIS/bioclim2070/hd45bi7018.tif")
bioclim2070_19 <- raster("./GIS/bioclim2070/hd45bi7019.tif")

bioclim2070 <- stack(bioclim2070_8,bioclim2070_9,bioclim2070_10,bioclim2070_11,bioclim2070_18,bioclim2070_19)
Variable Description
BIO1 Annual Mean Temperature
BIO2 Mean Diurnal Range (Mean of monthly (max temp – min temp))
BIO3 Isothermality (BIO2/BIO7) (x 100)
BIO4 Temperature Seasonality (standard deviation x 100)
BIO5 Max Temperature of Warmest Month
BIO6 Min Temperature of Coldest Month
BIO7 Temperature Annual Range (BIO5-BIO6)
BIO8 Mean Temperature of Wettest Quarter
BIO9 Mean Temperature of Driest Quarter
BIO10 Mean Temperature of Warmest Quarter
BIO11 Mean Temperature of Coldest Quarter
BIO12 Annual Precipitation
BIO13 Precipitation of Wettest Month
BIO14 Precipitation of Driest Month
BIO15 Precipitation Seasonality (Coefficient of Variation)
BIO16 Precipitation of Wettest Quarter
BIO17 Precipitation of Driest Quarter
BIO18 Precipitation of Warmest Quarter
BIO19 Precipitation of Coldest Quarter

Crop climate data to study area extent:

alps_bioclim <- crop(bioclim, extent(alps_boundary)+1)
alps_bioclim2070 <- crop(bioclim2070, extent(alps_dem90))

extent(alps_bioclim)
extent(alps_bioclim2070)

Resample environmental data and save object:

alps_bioclim <- resample(alps_bioclim, slope_aggregate)
extent(alps_bioclim)

writeRaster(alps_bioclim,filename="./GIS/alps_bioclim.grd", overwrite=T)
alps_bioclim <- stack("./GIS/alps_bioclim.grd")
dim(alps_bioclim)

alps_lacov <- resample(alps_lacov, alps_bioclim)
writeRaster(alps_lacov,filename="./GIS/alps_lacov.grd", overwrite=T)
alps_lacov <- stack("./GIS/alps_lacov.grd")

alps_bioclim2070 <- resample(alps_bioclim2070, slope_aggregate)
writeRaster(alps_bioclim2070,filename="./GIS/alps_bioclim2070.grd", overwrite=T)
alps_bioclim2070 <- stack("./GIS/alps_bioclim2070.grd")

Create an object including all environmental variables:

env <- stack(alps_bioclim, slope_aggregate, aspect_aggregate, alps_lacov, alps_dist)
env2070 <- stack(alps_bioclim2070, slope_aggregate, aspect_aggregate, alps_lacov, alps_dist)

names(env)
env@layers[[9]]@data@names<-"landcover"
env@layers[[10]]@data@names<-"distance"

names(env2070)
env2070@layers[[1]]@data@names<-"bio8_16"
env2070@layers[[2]]@data@names<-"bio9_16"
env2070@layers[[3]]@data@names<-"bio10_16"
env2070@layers[[4]]@data@names<-"bio11_16"
env2070@layers[[5]]@data@names<-"bio18_16"
env2070@layers[[6]]@data@names<-"bio19_16"
env2070@layers[[9]]@data@names<-"landcover"
env2070@layers[[10]]@data@names<-"distance"

env_mask <- raster::mask(env, mask=alps_ulrr_buf)
env2070_mask <- raster::mask(env2070, mask=alps_ulrr_buf)

writeRaster(env,filename="./GIS/env.grd", overwrite=T)
writeRaster(env_mask,filename="./GIS/env_mask.grd", overwrite=T)
writeRaster(env2070,filename="./GIS/env2070.grd", overwrite=T)
writeRaster(env2070_mask,filename="./GIS/env2070_mask.grd", overwrite=T)

Load environmental data:

env <- stack("./GIS/env.grd")
env2070 <- stack("./GIS/env2070.grd")
env_mask <- stack("./GIS/env_mask.grd")
env2070_mask <- stack("./GIS/env2070_mask.grd")

Load species data. Generate points out of polygons (according to Aulagnier et al. (2014) 31670 individuals are found in the Alps):

ibex <- readOGR("./GIS/Capraibex_iucn/species_42397.shp", layer="species_42397")
ibex_points <- spsample(ibex, n = 31670, "random")
ibex_sp <- SpatialPoints(coords=ibex_points, proj4string=CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

Visualize data: Digital elevation model and species polygons.

plot(alps_boundary)
plot(ibex, add=T)

plot of chunk unnamed-chunk-19

Visualize data: Digital elevation model and randomly computed species points.

plot(alps_boundary)
plot(ibex_sp, cex=0.01, add=T)

plot of chunk unnamed-chunk-20

Select species occurrence points where environmental information is available:

species_mask <- ibex_sp[complete.cases(extract(env_mask, ibex_sp)), ]

Test for collinearity and visualization:

cotest <- cor(getValues(env_mask), use = "complete.obs")
plotcorr(cotest, col=ifelse(abs(cotest) > 0.7,"red","grey"))

Subset of nearly uncorrelated environmental variables, if necessary:

env_mask <- subset(env_mask, c("bio10_16", "bio19_16","slope","aspect", "landcover", "distance"))
## Error: 'subset' must be logical
env2070_mask <- subset(env2070_mask, c("bio10_16", "bio19_16","slope","aspect", "landcover", "distance"))
## Error: 'subset' must be logical
cotest <- cor(getValues(env_mask), use = "complete.obs")
plotcorr(cotest, col=ifelse(abs(cotest) > 0.7,"red","grey"))

plot of chunk unnamed-chunk-24

Check data for NAs.

summary(env_mask)
na.omit(env_mask@data)

Selecting presence and background values:

set.seed(2)
background <- randomPoints(env_mask, 2000, species_mask)
presence <- gridSample(species_mask, env_mask, n = 1)

Combining presence and background points:

fulldata <- SpatialPointsDataFrame(rbind(presence, background),data = data.frame("species" = rep(c(1,0), c(nrow(presence), nrow(background)))),match.ID = FALSE,proj4string = CRS(projection(env_mask)))

Add information of environmental conditions at point locations:

fulldata@data <- cbind(fulldata@data, extract(env_mask, fulldata))

Check data for NAs.

Split data set into a training and test data set:

set.seed(2)
fold <- kfold(fulldata, k = 5)
traindata <- fulldata[fold != 1, ]
testdata <- fulldata[fold == 1, ]

Vector including variable names:

varnames <- c("bio10_16","bio19_16", "slope","aspect", "landcover", "distance")

Modelling

Generalized Additive Model:
gammodel <- gam(species ~ s(bio10_16) + s(bio19_16) + s(slope) + s(aspect) + s(landcover) + s(distance), family="binomial", data=traindata) 

summary(gammodel) 
plot(gammodel)

Model evaluation on test data:

a) Predict to test data

gamtest <- predict(gammodel, newdata = testdata, type = "response")

b) Calculate performance indices

val.prob(gamtest, testdata[["species"]])

plot of chunk unnamed-chunk-34

##        Dxy    C (ROC)         R2          D   D:Chi-sq        D:p 
##  8.940e-01  9.470e-01  6.858e-01  3.493e-01  1.778e+03         NA 
##          U   U:Chi-sq        U:p          Q      Brier  Intercept 
## -3.930e-04 -6.821e-13  1.000e+00  3.497e-01  2.648e-02 -1.523e-01 
##      Slope       Emax        S:z        S:p       Eavg 
##  1.092e+00  4.624e-02 -2.189e+00  2.857e-02  6.190e-03

Variable importance:

gamimp <- varImpBiomod(gammodel, varnames, na.omit(traindata@data))
barplot(100*gamimp/sum(gamimp), ylab = "Variable importance (%)", ylim=c(0,100))

plot of chunk unnamed-chunk-35

Prediction map. Probability of species occurrence for current and future data (2070, rcp45, HadGEM2-AO) :

gammap <- predict(env_mask, gammodel, type = "response")
plot(gammap)
plot(alps_boundary, add=T)
plot(ibex_sp,cex=0.01, add=T)
title("Actual")
mtext("Longitude [°]",side=1,line=-3,cex=0.8,outer=TRUE)
mtext("Latitude [°]",side=2,line=-1.3,cex=0.8,outer=TRUE)
north.arrow(yb=43.5, xb=16, len=0.3,lab="N",cex.lab=1, col="black")

plot of chunk unnamed-chunk-36

gammap2070 <- predict(env2070_mask, gammodel, type = "response")
plot(gammap2070)
plot(alps_boundary, add=T)
plot(ibex_sp,cex=0.01, add=T)
title("Future")
mtext("Longitude [°]",side=1,line=-3,cex=0.8,outer=TRUE)
mtext("Latitude [°]",side=2,line=-1.3,cex=0.8,outer=TRUE)
north.arrow(yb=43.5, xb=16, len=0.3,lab="N",cex.lab=1, col="black")

plot of chunk unnamed-chunk-36

Random Forest:
rftraindata <- as(traindata, "data.frame")
rftraindata$species <- as.factor(rftraindata$species)
rftraindata <- rftraindata[complete.cases(rftraindata), ]

rfmodel <- randomForest(species ~  bio10_16 + bio19_16 + slope + aspect + landcover + distance , data = rftraindata)

Model evaluation on test data:

a) Predict to test data

rftest <- predict(rfmodel, newdata = testdata, type = "prob")[,2]

b) Calculate performance indices

val.prob(rftest, testdata[["species"]])

plot of chunk unnamed-chunk-39

##        Dxy    C (ROC)         R2          D   D:Chi-sq        D:p 
##  9.053e-01  9.527e-01  7.330e-01  3.954e-01  1.650e+03         NA 
##          U   U:Chi-sq        U:p          Q      Brier  Intercept 
## -4.795e-04 -1.819e-12  1.000e+00  3.959e-01  1.984e-02 -2.729e-01 
##      Slope       Emax        S:z        S:p       Eavg 
##  1.230e+00  8.861e-02 -3.851e+00  1.177e-04  1.303e-02

Variable importance:

rfImp <- importance(rfmodel)
varImpPlot(rfmodel)

plot of chunk unnamed-chunk-40

Response functions:

par(mfrow=c(2,3))
for (i in seq_along(varnames)) {partialPlot(rfmodel, rftraindata, varnames[i], xlab = varnames[i], main="")}
## Error: object 'varnames' not found

Prediction map:

rfmap <- predict(env_mask, rfmodel, type = "prob", na.rm=F, index=2)
plot(rfmap)
plot(alps_boundary, add=T)
plot(ibex_sp,cex=0.01, add=T)
title("Actual")
mtext("Longitude [°]",side=1,line=-3,cex=0.8,outer=TRUE)
mtext("Latitude [°]",side=2,line=-1.3,cex=0.8,outer=TRUE)
north.arrow(yb=43.5, xb=16, len=0.3,lab="N",cex.lab=1, col="black")