#################################################### ###### Species Distribution Modeling tutorial ###### #################################################### # created for LifeWatch, November 2014 # by Mirko Di Febbraro (mirkodifebbraro@gmail.com) # Departement of Biosciences and Territory (DiBT), University of Molise. Italy # and # by Paolo Colangelo (paolo.colangelo@uniroma1.it) # Institute of Ecosistem Study (ISE), National Research Council, Italy ##### MODULE 1 ##### # load the library dismo library(dismo) #download the occurnces of P. krameri krameri <- gbif(genus = "Psittacula", species = "krameri", geo = T, concept = T) # check variable names colnames(krameri) #We want only lat, long and coordinates Uncertainty krameri <- krameri[, c(7:9)] #check the new dataset head(krameri) str(krameri) summary(krameri) # convert coordinates uncertainty column into numeric vector krameri$coordUncertaintyM <- as.numeric(krameri$coordUncertaintyM) #exclude individuals with NA's in coordinates krameri <- krameri[!is.na(krameri$lat), ] # remove duplicated records with identical coordinates xy <- cbind(krameri$lon, krameri$lat) dup <- duplicated(xy) krameri <- krameri[!dup, ] #check dataset consistency dim(krameri) # plot the records library(maptools) data(wrld_simpl) plot(wrld_simpl) points(krameri$lon, krameri$lat, pch = 21, bg = "red", cex = 0.5) # remove occurrences falling outside the range of the studied species. Import range shapefile with the readShapePoly function (maptools package) krameri_EOO <- readShapePoly("Shapefile_Pkrameri/Psittacula_krameri_1529_BL.shp") krameri_EOO <- krameri_EOO[krameri_EOO@data$ORIGIN == 1, ] plot(wrld_simpl) # plot species range plot(krameri_EOO, col = "cyan", add = T) xy <- cbind(krameri$lon, krameri$lat) points(xy, pch = 21, bg = "red", cex = 0.5) # remove occurrences not falling into the species range # First convert krameri into a sp object krameri_sp <- SpatialPoints(xy) # overlay krameri occurrences to species range ov <- over(krameri_sp, krameri_EOO) # remove rows with NAs (not overlaying to species range) krameri <- krameri[!is.na(ov[, 1]), ] plot(wrld_simpl) plot(krameri_EOO, col = "cyan", add = T) xy <- cbind(krameri$lon, krameri$lat) points(xy, pch = 21, bg = "red", cex = 0.5) dim(krameri) write.table(krameri, file = "krameri_NativeRange.txt") ###################### #### END MODULE 1 #### ###################### #################### ##### MODULE 2 ##### #################### # reload the file obtained in module 1 krameri <- read.table(file = "krameri_NativeRange.txt", header = T) library(raster) # download bioclim variables from worldclim dataset at a spatial resolution of 2.5 minutes of a degree r <- getData("worldclim", var = "bio", res = 2.5) #plot(r$bio1) ##### Remove duplicates occurences within the same cell ##### occ.desaggregation_RASTER <- function(df, colxy, colvar, rast, plot = T) { require(raster) df_ini <- df var1 <- cellFromXY(rast, df[, colxy]) if (any(is.na(var1))) { stop("no NA admitted in species occurrences!") } var2 <- split(var1, var1) l <- sapply(var2, length) l1 <- l[l > 1] for (j in 1:length(l1)) { w <- which(var1 %in% as.numeric(names(l1[j]))) found <- df[w, ] s <- sample(w, 1) df[w[w != s], colxy] <- NA } df_final <- df[!is.na(df[, colxy[1]]), ] if (plot == T) { x11() plot(df_ini[, colxy], main = "distribution of occurences", sub = paste("# initial (black):", nrow(df_ini), " | # kept (red): ", nrow(df_final)), pch = 19, col = "black", cex = 0.5) points(df_final[, colxy], pch = 19, col = "red", cex = 0.2) } return(df_final) } krameri <- occ.desaggregation_RASTER(df = krameri, colxy = c(2, 1), rast = r[[1]], plot = F) dim(krameri) ##### Reducing the sampling bias with subsampling procedures ##### #visual check for spatial bias library(maptools) data(wrld_simpl) plot(wrld_simpl) points(krameri$lon, krameri$lat, pch = 21, bg = "red", cex = 0.5) # subsample the occurrences (function provided with the tutorial) # for a description of the function's arguments see head of script source("script/spatial.bias.reduction.R") reduction <- spatial.bias.reduction(points = krameri, xy.cols = c(2, 1), grid.cell.size = 100, mask_base = r[[1]], FUN = "median", sel.by.err = F, rep_bsp = 1) # plot species range and countries boundaries krameri_EOO <- readShapePoly("Shapefile_Pkrameri/Psittacula_krameri_1529_BL.shp") krameri_EOO <- krameri_EOO[krameri_EOO@data$ORIGIN == 1, ] plot(krameri_EOO) plot(wrld_simpl, add = T) plot(krameri_EOO, add = T, col = "cyan") # plot grid used for the subsampling procedure plot(reduction[[2]], add = T) # plot occurrences pre (blue) and post (red) subsampling procedure points(krameri[c(2, 1)], pch = 19, col = "blue", cex = 0.5) points(reduction[[1]][c(2, 1)], pch = 19, col = "red", cex = 0.5) krameri <- reduction[[1]] dim(krameri) ###### Coordinates uncertainty: how to deal with this issue ####### # check for smissing value summary(krameri$coordUncertaintyM) # choose 100 km as maximum error for records without information on the spatial uncertainty krameri[is.na(krameri[, 3]), 3] <- 100000 source("script/coordinates.precision.resampling.R") resampling <- coordinates.precision.resampling(species = krameri, xy.cols = c(2, 1), err = 2500, err.col = 3, vars = r[[1]], rep_cp = 3) plot(krameri_EOO, col = "cyan") plot(wrld_simpl, add = T, border = "grey", col = "grey") plot(krameri_EOO, col = "cyan", add = T, border = "grey") plot(SpatialPoints(krameri[c(2, 1)]), pch = 17, add = T, cex = 0.7) plot(SpatialPoints(resampling[[1]]), add = T, col = "red") plot(SpatialPoints(resampling[[2]]), add = T, col = "blue") plot(SpatialPoints(resampling[[3]]), add = T, col = "darkgreen") krameri <- resampling #the new dataset included 3 replicartes of 197 occurences each. This is the dataset that we will use to train the model in section 3. Save it! save(krameri, file="krameri.training") ### because kramery is a list with 3 replicates we will save it as R object. You can open for section 3 using the function load() ### create a 1000 km buffer around occurrences to be used as model training area library(rgeos) buff<-gBuffer(SpatialPoints(do.call(rbind, krameri)), width=res(r)[1]*100) # resolution of r is almost 5 km --> 5km*100=500km vars<-crop(r, buff, progress="text") vars<-mask(vars, buff, progress="text") plot(vars[[1]]) plot(wrld_simpl, add = T, border = "grey", col = "grey") plot(vars[[1]], add=T) ##### Dealing with multi-collinearity: Pearson's r and Variance Inflation Factor (VIF) analysis ##### ### Pearson's r ### # extract values of environmental raster at 10000 random cells sample<-randomPoints(vars, 10000) vars_tab<-extract(vars, sample) # check which variables have a Pearson's r greater than 0.7 abs(cor(vars_tab)) > 0.7 # select predictors to be removed abs(cor(vars_tab[-c(3,4,5,6,8,9,10,11,13,14,16,18), -c(3,4,5,6,8,9,10,11,13,14,16,18)])) > 0.7 # create the new set of rasters vars_r <- vars[[-c(3,4,5,6,8,9,10,11,13,14,16,18)]] ### variance inflation factor analisys (function provided with the tutorial) ### source("script/vif_functions.R") # the function find and remove each predictor reporting a vif greater than 3 vars_vif <- vif(vars_tab, thresh = 3, trace = F) # create the new set of rasters vars_vif <- vars[[which(names(vars) %in% vars_vif)]] # there are some differences in the results names(vars_r) names(vars_vif) dir.create("predictors") writeRaster(vars_r, "predictors/", bylayer=T, suffix="names") #################### ##### MODULE 3 ##### #################### #load libraries library(rgeos) library(raster) library(maptools) library(gtools) library(biomod2) data(wrld_simpl) # load training dataset and training variables load("krameri.training") vars_training<-stack(mixedsort(list.files("predictors/", "gri", full.names=T))) # create raster layers for the model projection under the current climate condition at a global scale r <- getData("worldclim", var = "bio", res = 2.5) r<-r[[which(names(r)%in%names(vars_training))]] #downscale to 5Km vars_proj<-stack(aggregate(r, 5, progress="text")) # create raster layers for the model projection under future climate conditions at a global scale # download file download.file("http://biogeo.ucdavis.edu/data/climate/cmip5/2_5m/cc85bi70.zip", "future.zip") unzip("future.zip", exdir="future") #load layers future<-stack(mixedsort(list.files("future","tif",full.names=T))) names(future)<-paste("bio", 1:19, sep="") future<-future[[which(names(future)%in%names(vars_training))]] #downscale to 5Km future<-stack(aggregate(future, 5, progress="text")) #create output directiry a set the new working directory dir.create("biomod.output") setwd("biomod.output") # run the model with the 1st replicated dataset i=1 # set the input data for the model training mydata<-BIOMOD_FormatingData(resp.var=rep(1, nrow(krameri[[i]])), expl.var=vars_training, PA.nb.rep=1, PA.nb.absences=1000, PA.strategy='random', resp.xy=krameri[[i]], resp.name=paste("krameri.", i, sep="")) mydata # set the modelling options for the algorithms intended to be used myBiomodOption<-BIOMOD_ModelingOptions(GLM=list(type='polynomial'), GBM=list(n.trees=1000), GAM=list(k=4), MAXENT=list(path_to_maxent.jar="../", maximumiterations=1000)) myBiomodOption # set parameters for the model training mymodel<-BIOMOD_Modeling(data=mydata, models=c('GLM','GAM','GBM','MAXENT'), models.options=myBiomodOption, NbRunEval=3, DataSplit=70, VarImport=2, models.eval.meth=c('ROC', 'TSS'), do.full.models=F) mymodel #print out variable importance get_variables_importance(mymodel) # project models over current climate in the trainign area (it may takes a while) #downscale bioclimatic variables to save time vars_training_5km<-stack(aggregate(vars_training, 5, progress="text")) myproj<-BIOMOD_Projection(modeling.output=mymodel, new.env= vars_training_5km, proj.name="current_training", selected.models= "all", do.stack=T) # plot output probability maps of the 1st replicate (4 model, only the first run = 4 maps) source("../script/palette_MAXENT.R") map.current<-stack("krameri.1/proj_current_training/proj_current_training_krameri.1.gri") plot(map.current[[c(1,2,3,4)]], col=palette_MAXENT(min=0, max=1000, interval=100)[[1]], breaks=palette_MAXENT(0, 1000, 100)[[2]]) # project models under current climate conditions at the global scale (it may takes a while) myproj<-BIOMOD_Projection(modeling.output=mymodel, new.env= vars_proj, proj.name="current", selected.models= "all", do.stack=T) # plot output probability maps of the 1st replicate (4 model x 3 runs = 12 maps) source("../script/palette_MAXENT.R") map.current<-stack("krameri.1/proj_current/proj_current_krameri.1.gri") plot(map.current, col=palette_MAXENT(min=0, max=1000, interval=100)[[1]], breaks=palette_MAXENT(0, 1000, 100)[[2]]) # Create the ensemble model myEM<-BIOMOD_EnsembleModeling(modeling.output=mymodel, em.by="all", eval.metric=c('ROC','TSS'), eval.metric.quality.threshold=c(0.7,0.4), models.eval.meth = c('ROC', 'TSS'), prob.mean=F, prob.mean.weight=T) # create the probability and binary ensemble maps myEF<-BIOMOD_EnsembleForecasting(EM.output=myEM, projection.output=myproj, total.consensus=T, binary.meth=c("TSS")) # project models for future climate condition myproj<-BIOMOD_Projection(modeling.output=mymodel, new.env=future, proj.name="future", selected.models="all", do.stack=T) # create the probability and binary ensemble maps myEF<-BIOMOD_EnsembleForecasting(EM.output=myEM, projection.output=myproj, total.consensus=T, binary.meth=c("TSS")) # plot probability maps of the Ensemble models for the 1st replicate (current and future climatic conditions) source("../script/palette_MAXENT.R") map.current<-stack("krameri.1/proj_current/proj_current_krameri.1_ensemble.gri") plot(map.current[[2]], col=palette_MAXENT(min=0, max=1000, interval=100)[[1]], breaks=palette_MAXENT(0, 1000, 100)[[2]]) map.future<-stack("krameri.1/proj_future/proj_future_krameri.1_ensemble.gri") plot(map.future[[2]], col=palette_MAXENT(0, 1000, 100)[[1]], breaks=palette_MAXENT(0, 1000, 100)[[2]]) ### We build other ensemble models based on second and third re replicate for(i in 2:length(krameri)) { mydata<-BIOMOD_FormatingData(resp.var=rep(1, nrow(krameri[[i]])), expl.var=vars_training, PA.nb.rep=1, PA.nb.absences=1000, PA.strategy='random', resp.xy=krameri[[i]], resp.name=paste("krameri.", i, sep="")) mymodel<-BIOMOD_Modeling(data=mydata, models=c('GLM','GAM','GBM','MAXENT'), models.options=myBiomodOption, NbRunEval=3, DataSplit=70, VarImport=2, models.eval.meth=c('ROC', 'TSS'), do.full.models=F) myEM<-BIOMOD_EnsembleModeling(modeling.output=mymodel, em.by="all", eval.metric=c('ROC','TSS'), eval.metric.quality.threshold=c(0.7,0.4), models.eval.meth = c('ROC', 'TSS'), prob.mean=F, prob.mean.weight=T) myproj<-BIOMOD_Projection(modeling.output=mymodel, new.env=vars_proj, proj.name="current", selected.models="all", do.stack=T) myEF<-BIOMOD_EnsembleForecasting(EM.output=myEM, projection.output=myproj, total.consensus=T, binary.meth=c("TSS")) myproj<-BIOMOD_Projection(modeling.output=mymodel, new.env=future, proj.name="future", selected.models="all", do.stack=T) myEF<-BIOMOD_EnsembleForecasting(EM.output=myEM, projection.output=myproj, total.consensus=T, binary.meth=c("TSS")) } ### END MODULE 3 #################### ##### MODULE 4 ##### #################### ## load ensemble models files l<-list.files(getwd(), "ensemble.models.out", recursive=T, full.names=T) for(i in 1:length(l)) { load(l[i]) } myEM<-lapply(grep("ensemble.models.out", ls(), value=T), get) # get TSS values for the ensemble models library(biomod2) TSS<-unlist(lapply(myEM, get_evaluations), recursive=F) TSS<-sapply(TSS[grep("TSS", names(TSS))], "[[", c(2)) TSS # load ensemble probability maps m<-list.files(getwd(), "ensemble.gri", recursive=T, full.names=T) m_current<-lapply(grep("current", m, value=T), stack) m_future<-lapply(grep("future", m, value=T), stack) # and take ensembles obtained averaging the single models by TSS m_current<-lapply(m_current, "[[", c(2)) m_future<-lapply(m_future, "[[", c(2)) # compute the weigthed mean of the ensembles current_final<-weighted.mean(stack(m_current), w=TSS) future_final<-weighted.mean(stack(m_future), w=TSS) # plot the resulting maps source("script/palette_MAXENT.R") plot(current_final, col=palette_MAXENT(min=0, max=1000, interval=100)[[1]], breaks=palette_MAXENT(0, 1000, 100)[[2]]) plot(future_final, col=palette_MAXENT(min=0, max=1000, interval=100)[[1]], breaks=palette_MAXENT(0, 1000, 100)[[2]]) # load binary maps b<-list.files(getwd(), "ensemble_TSSbin.gri", recursive=T, full.names=T) b_current<-lapply(grep("current", b, value=T), stack) b_future<-lapply(grep("future", b, value=T), stack) # retain just ensembles made by using TSS b_current<-lapply(b_current, "[[", c(2)) b_future<-lapply(b_future, "[[", c(2)) # consider only cells of predicted presence by at least 2 of the 3 replicates and make the sum b_current<-sum(stack(b_current)) values(b_current)[values(b_current)%in%c(0:1)]<-0 values(b_current)[values(b_current)%in%c(2:3)]<-1 b_future<-sum(stack(b_future)) values(b_future)[values(b_future)%in%c(0:1)]<-0 values(b_future)[values(b_future)%in%c(2:3)]<-1 # plot the resulting maps plot(b_current, col=c("darkgrey", "darkgreen")) plot(b_future, col=c("darkgrey", "darkgreen")) # evaluate the effect of climate change library(scales) par(mfrow=c(1,3)) hist(values(current_final), xlim=c(0,1000), xlab="Probability of presence*1000", col="darkblue", main="Predicted probabilies of presence") hist(values(future_final), xlim=c(0,1000), xlab="Probability of presence*1000", col=alpha("red",0.7), main="Future climate", add=T) barplot(c(table(values(b_current))[2], table(values(b_future))[2]), col="darkgreen", main="Area of predicted presence", names=c("Current", "Future"), ylab="number of cells") rtp_current<-rasterToPoints(b_current) rtp_current<-rtp_current[rtp_current[,3]==1,2] rtp_future<-rasterToPoints(b_future) rtp_future<-rtp_future[rtp_future[,3]==1,2] boxplot(list(rtp_current, rtp_future), col="grey", main="Latitude of predicted presence", names=c("Current", "Future"), ylab="decimal degrees") # evaluate the effect of extrapolation in model projections (just for the 1st replicate) ## load single models files s<-list.files(getwd(), "formated.input.data", recursive=T, full.names=T, all.files=T) load(s[1]) s<-grep("data", ls(), value=T) mydata<-get(s) # extract training dataset dataset<-mydata@coord # perform MESS analysis (it may takes several minutes) library(gtools) library(dismo) r <- getData("worldclim", var = "bio", res = 2.5) vars_training<-stack(mixedsort(list.files("predictors/", "gri", full.names=T))) r<-r[[which(names(r)%in%names(vars_training))]] vars_proj<-stack(aggregate(r, 5, progress="text")) MESS<-mess(x=vars_proj, v=extract(vars_training, dataset), full=F) color.high<-colorRampPalette(c("lightgrey","blue")) color.low<-colorRampPalette(c("red","lightgrey")) plot(MESS, col=c(color.low(abs(minValue(MESS))),color.high(maxValue(MESS))))