library(biomod2) set.seed(22) # load species data DataSpecies <- read.csv(file="Current Tufted Puffin Distribution.csv", header=TRUE) head(DataSpecies) # the name of studied species myRespName <- 'Tufted.Puffin' # the presence/absences data for our species myResp <- as.numeric(DataSpecies[,myRespName]) # the XY coordinates of species data myRespXY <- DataSpecies[,c("LONGITUDE","LATITUDE")] # Environmental variables extracted from Worldclim myExpl<-stack(getData('worldclim', var='bio', res=5)) myExpl<-stack(crop(myExpl,(extent(-180, -120, 33, 69)))) names(myExpl)<-c("Annual Mean Temp", "Mean Diurnal Range", "Isothermality", "Temp Seasonality", "Max Temp of Warmest", "Min Temp of Coldest", "Temp Annual Range", "Mean Temp of Wettest Quarter", "Mean Temp of Driest Quarter", "Mean Temp of Warmest Quarter", "Mean Temp of Coldest Quarter", "Annual Precip", "Precip of Wettest Month", "Precip of Driest Month", "Precip Seasonality", "Precip of Wettest Quarter", "Precip of Driest Quarter", "Precip of Warmest Quarter", "Precip of Coldest Quarter") #creation of distance to ocean layer (helps projections restrict to shoreline) z<-raster(xmn= -180, xmx= -120, ymn= 33, ymx= 69, nrow= 432, ncol= 720) z[]= runif(432*720) zz<-mask(z, myExpl$Mean.Temp.of.Warmest.Quarter, inverse=T) dist.layer<-distance(zz, doEdge=T) dist.layer<-dist.layer/1000 s <- calc(dist.layer, fun=function(x){ x[x > 200] <- NA; return(x)} ) myExpl<-stack((myExpl$Temp.Annual.Range/10), (myExpl$Mean.Diurnal.Range/10), (myExpl$Mean.Temp.of.Warmest.Quarter/10), (myExpl$Annual.Precip/100), (myExpl$Precip.of.Warmest.Quarter/100), (dist.layer)) myExpl<-mask(myExpl, s) myExpl<-stack(myExpl) names(myExpl)<-c("ATR", "MDR", "MTWQ", "AP", "PWQ", "DIST") plot(myExpl) ### formating_data myBiomodData <- BIOMOD_FormatingData(resp.var = myResp, expl.var = myExpl, resp.xy = myRespXY, resp.name = myRespName, PA.nb.rep = 2, PA.nb.absences = 1000, PA.strategy = "disk", PA.dist.min = 30000, PA.dist.max = NULL) ### print_formating_data myBiomodData ### plot_formating_data plot(myBiomodData) level.plot(data.in=myResp, XY= myRespXY) ### modeling_options #Defining Models Options using default options myBiomodOption <- BIOMOD_ModelingOptions( GLM = list( type = 'quadratic', interaction.level = 0, myFormula = NULL, test = 'AIC', family = binomial(link = 'logit'), mustart = 0.5, control = glm.control(epsilon = 1e-08, maxit = 50, trace = FALSE) ), GBM = list( distribution = 'bernoulli', n.trees = 2500, interaction.depth = 7, n.minobsinnode = 5, shrinkage = 0.001, bag.fraction = 0.5, train.fraction = 1, cv.folds = 3, keep.data = FALSE, verbose = FALSE, perf.method = 'cv'), RF = list( do.classif = TRUE, ntree = 500, mtry = 'default', nodesize = 5, maxnodes = NULL) ) ### modeling myBiomodModelOut <- BIOMOD_Modeling(data= myBiomodData, models = c("GLM", "GBM", "RF"), models.options = myBiomodOption, NbRunEval= 20, DataSplit= 70, Prevalence=0.5, VarImport=2, models.eval.meth = c('TSS','ROC', 'ACCURACY'), SaveObj = TRUE, rescal.all.models = FALSE, do.full.models = FALSE, modeling.id=paste(myRespName, "First Modeling", sep="")) ### modeling_summary myBiomodModelOut ### modeling model evaluation # get all models evaluation myBiomodModelEval <- get_evaluations(myBiomodModelOut, as.data.frame=T) myBiomodModelEval GLM.eval<-myBiomodModelEval[grep("GLM", myBiomodModelEval$Model.name),] TSS.GLM.eval<-GLM.eval[grep("TSS", GLM.eval$Eval.metric),] AUC.GLM.eval<-GLM.eval[grep("ROC", GLM.eval$Eval.metric),] TSS.GLM.mean<-mean(TSS.GLM.eval$Testing.data) TSS.GLM.sd<-sd(TSS.GLM.eval$Testing.data) AUC.GLM.mean<-mean(AUC.GLM.eval$Testing.data) AUC.GLM.sd<-sd(AUC.GLM.eval$Testing.data) GBM.eval<-myBiomodModelEval[grep("GBM", myBiomodModelEval$Model.name),] TSS.GBM.eval<-GBM.eval[grep("TSS", GBM.eval$Eval.metric),] AUC.GBM.eval<-GBM.eval[grep("ROC", GBM.eval$Eval.metric),] TSS.GBM.mean<-mean(TSS.GBM.eval$Testing.data) TSS.GBM.sd<-sd(TSS.GBM.eval$Testing.data) AUC.GBM.mean<-mean(AUC.GBM.eval$Testing.data) AUC.GBM.sd<-sd(AUC.GBM.eval$Testing.data) RF.eval<-myBiomodModelEval[grep("RF", myBiomodModelEval$Model.name),] TSS.RF.eval<-RF.eval[grep("TSS", RF.eval$Eval.metric),] AUC.RF.eval<-RF.eval[grep("ROC", RF.eval$Eval.metric),] TSS.RF.mean<-mean(TSS.RF.eval$Testing.data) TSS.RF.sd<-sd(TSS.RF.eval$Testing.data) AUC.RF.mean<-mean(AUC.RF.eval$Testing.data) AUC.RF.sd<-sd(AUC.RF.eval$Testing.data) # print the dimnames of this object dimnames(myBiomodModelEval) ### modeling variable importance # print variable importances var.imp<-get_variables_importance(myBiomodModelOut, as.data.frame=T) var.imp data.temp <- cbind(TOTAL=get_formal_data(myBiomodModelOut,'resp.var'), get_formal_data(myBiomodModelOut,'expl.var')) library(plyr) Pres.abs<-factor(data.temp[,1]) Pres.abs<-revalue(Pres.abs, c("0"="Absence", "1"="Presence")) pca <- princomp(data.temp[,-1], cor= T) # principal components analysis using correlation matrix loadings(pca) summary(pca) pc.comp <- pca$scores pc.comp<-pca$scores[1:781,] loadings<-loadings(pca) loadings Comp.1<-as.numeric(loadings[,1]) Comp.1<-round(Comp.1, digits=3) Comp.2<-as.numeric(loadings[,2]) Comp.2<-round(Comp.2, digits=3) PCAs<-cbind(Comp.1, Comp.2) abcd<-data.frame(PCAs) row.names(abcd)<-c("ATR", "MDR", "MTWQ", "AP", "PWQ", "DIST") # k-means clustering [assume 3 clusters] km <- kmeans(pc.comp, centers=2, nstart=5) ggdata <- data.frame(pc.comp, Cluster=km$cluster, Species=Pres.abs[1:781]) png(filename="Figure 2", width=1200, height=900) #PCA plot library(gridBase) library(gridExtra) tt.2<-ttheme_minimal(base_size=11) ggplot(ggdata) + geom_point(aes(x=Comp.1, y=Comp.2, color=factor(Species)), size=3, shape=20) + stat_ellipse(aes(x=Comp.1,y=Comp.2,fill=factor(Species)), geom="polygon", level=0.95, alpha=0.2) + guides(color=guide_legend("Species"), fill=guide_legend("Species"))+ theme(legend.title=element_blank()) + #theme(legend.position = c(.94, .5))+ theme(legend.background = element_rect(fill=alpha("Blue", 0.0 )))+ theme(legend.text =element_text(size=16))+ #ggtitle("Principal Component Analysis")+ theme(plot.title = element_text(size = 14, lineheight=1, face="bold"))+ annotation_custom(tableGrob(d=abcd, theme=tt.2), xmin=4.5, xmax=5.7, ymin=-3, ymax=-1.5) #ggExtra::ggMarginal(pca.plot, data= ggdata, type = "density", margins = "y", size = 7) dev.off() ### ensemble modeling myBiomodEM <- BIOMOD_EnsembleModeling(modeling.output = myBiomodModelOut, chosen.models = 'all', em.by='all', eval.metric = c('TSS'), eval.metric.quality.threshold = c(0.7), prob.mean = F, prob.cv = F, prob.ci = F, prob.ci.alpha = 0.05, prob.median = F, committee.averaging = F, prob.mean.weight = T, prob.mean.weight.decay = 'proportional' ) ### ensemble_modeling outputs # print summary myBiomodEM list.files(myBiomodEM) # get evaluation scores get_evaluations(myBiomodEM) plot(myBiomodEM) ### projection curent # projection over the globe under current conditions myBiomodProj <- BIOMOD_Projection( modeling.output = myBiomodModelOut, new.env = myExpl, proj.name = 'present', selected.models = 'all', binary.meth = 'TSS', filtered.meth = "TSS", compress = 'xy', clamping.mask = F, output.format = '.grd') # summary of crated oject myBiomodProj # files created on hard drive list.files(myBiomodProj) ### projection curent plot # make some plots sub-selected by str.grep argument plot(myBiomodProj, str.grep = 'GLM') plot(myBiomodProj, str.grep = 'RF') plot(myBiomodProj, str.grep= "GBM") ### projection curent getProj # if you want to make custom plots, you can also get the projected map myCurrentProj <- get_predictions(myBiomodProj) myCurrentProj myCurrentProj.1<-mean(myCurrentProj/1000) myCurrentProj.1 ### Response curve analysis # # 4.1 Load the models for which we want to extract the predicted response curves myGLMs <- BIOMOD_LoadModels(myBiomodModelOut, models=c('GLM')) myRFs<-BIOMOD_LoadModels(myBiomodModelOut, models=c('RF')) myGBMs<-BIOMOD_LoadModels(myBiomodModelOut, models="GBM") #4.2 plot 2D response plots myRespPlot2D <- response.plot2(models = myGLMs, Data = get_formal_data(myBiomodModelOut, 'expl.var'), show.variables= get_formal_data(myBiomodModelOut, 'expl.var.names'), do.bivariate = F, display_title= F, fixed.var.metric = 'mean', legend = F, data_species = get_formal_data(myBiomodModelOut, 'resp.var')) library(dplyr) library(tidyr) resp.plot.dat <- bind_rows(lapply(myRespPlot2D, function(dat_){ dat_$id <- rownames(dat_) expl.dat_ <- dat_ %>% select(1, id) %>% gather("expl.name", "expl.val", 1) pred.dat_ <- dat_ %>% select(-1, id) %>% gather("pred.name", "pred.val", -id) out.dat_ <- full_join(expl.dat_, pred.dat_) return(out.dat_) })) head(resp.plot.dat) resp.ap<-resp.plot.dat gg.AP <- ggplot(resp.plot.dat[resp.plot.dat$expl.name=="AP",], aes(x = expl.val, y = pred.val, col = pred.name)) + geom_line()+ theme(legend.position = "none") + ylab("Predicited Suitability Value") + theme(plot.title= element_blank(), axis.title.x= element_blank()) + annotate("text", x=42, y=.90, label="D", size= 15) gg.ATR <- ggplot(resp.plot.dat[resp.plot.dat$expl.name=="ATR",], aes(x = expl.val, y = pred.val, col = pred.name)) + geom_line()+ theme(legend.position = "none") + ylab("Predicited Suitability Value") + theme(plot.title= element_blank(), axis.title.x= element_blank()) + annotate("text", x=48, y=.90, label="A", size= 15) gg.MDR <- ggplot(resp.plot.dat[resp.plot.dat$expl.name=="MDR",], aes(x = expl.val, y = pred.val, col = pred.name)) + geom_line()+ theme(legend.position = "none") + theme(plot.title= element_blank(), axis.title.x= element_blank(), axis.title.y= element_blank()) + annotate("text", x=18, y=.82, label="B", size= 15) gg.MTWQ <- ggplot(resp.plot.dat[resp.plot.dat$expl.name=="MTWQ",], aes(x = expl.val, y = pred.val, col = pred.name)) + geom_line()+ theme(legend.position = "none") + theme(plot.title= element_blank(), axis.title.x= element_blank(), axis.title.y= element_blank()) + annotate("text", x=22, y=.82, label="C", size= 15) gg.PWQ <- ggplot(resp.plot.dat[resp.plot.dat$expl.name=="PWQ",], aes(x = expl.val, y = pred.val, col = pred.name)) + geom_line()+ theme(legend.position = "none") + theme(plot.title= element_blank(), axis.title.x= element_blank(), axis.title.y= element_blank()) + annotate("text", x=7.75, y=.80, label="E", size= 15) gg.DIST <- ggplot(resp.plot.dat[resp.plot.dat$expl.name=="DIST",], aes(x = expl.val, y = pred.val, col = pred.name)) + geom_line()+ theme(legend.position = "none") + theme(plot.title= element_blank(), axis.title.x= element_blank(), axis.title.y= element_blank()) + annotate("text", x=180, y=.90, label="F", size= 15) grid.arrange(arrangeGrob(gg.ATR, gg.MDR, gg.MTWQ, ncol=3), arrangeGrob(gg.AP, gg.PWQ, gg.DIST, ncol=3), nrow=2) ### future projection # load environmental variables for the future #2050 RCP 4.5 myExplFuture4.5.1<- getData('CMIP5', var='bio', res=5, rcp=45, model='GS', year=50) myExplFuture4.5.2<-getData('CMIP5', var='bio', res=5, rcp=45, model='HD', year=50) myExplFuture4.5.3<- getData('CMIP5', var='bio', res=5, rcp=45, model='GF', year=50) myExplFuture4.5.4<- getData('CMIP5', var='bio', res=5, rcp=45, model='IP', year=50) myExplFuture4.5.5<- getData('CMIP5', var='bio', res=5, rcp=45, model='MP', year=50) myExplFuture4.5.6<- getData('CMIP5', var='bio', res=5, rcp=45, model='NO', year=50) myExplFuture4.5.7<- getData('CMIP5', var='bio', res=5, rcp=45, model='CC', year=50) myExplFuture4.5.8<- getData('CMIP5', var='bio', res=5, rcp=45, model='BC', year=50) myExplFuture4.5<-stack(myExplFuture4.5.1, myExplFuture4.5.2, myExplFuture4.5.3, myExplFuture4.5.4, myExplFuture4.5.5, myExplFuture4.5.6, myExplFuture4.5.7, myExplFuture4.5.8) myExplFuture4.5<-subset(myExplFuture4.5, c(1,20,39,58,77,96,115,134,2,21,40,59,78,97,116,135,3,22,41,60,79,98,117,136, 4,23,42,61,80,99,118,137,5,24,43,62,81,100,119,138,6,25,44,63,82,101,120,139,7,26,45,64,83,102,121,140,8,27,46,65,84,103,122,141, 9,28,47,66,85,104,123,142,10,29,48,67,86,105,124,143,11,30,49,68,87,106,125,144,12,31,50,69,88,107,126,145,13,32,51,70,89,108,127,146, 14,33,52,71,90,109,128,147,15,34,53,72,91,110,129,148,16,35,54,73,92,111,130,149,17,36,55,74,93,112,131,150,18,37,56,75,94,113,132,151, 19,38,57,76,95,114,133,152)) myExplFuture4.5<-stack(crop(myExplFuture4.5,(extent(-180, -120, 33, 69)))) myExplFuture4.5<-stackApply(myExplFuture4.5, c(1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5, 6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,8,8,8,8,8,8,8,8,9,9,9,9,9,9,9,9,10,10,10,10,10,10,10,10, 11,11,11,11,11,11,11,11,12,12,12,12,12,12,12,12,13,13,13,13,13,13,13,13,14,14,14,14,14,14,14,14, 15,15,15,15,15,15,15,15,16,16,16,16,16,16,16,16,17,17,17,17,17,17,17,17, 18,18,18,18,18,18,18,18,19,19,19,19,19,19,19,19), fun= mean) names(myExplFuture4.5)<-c("Annual Mean Temp", "Mean Diurnal Range", "Isothermality", "Temp Seasonality", "Max Temp of Warmest", "Min Temp of Coldest", "Temp Annual Range", "Mean Temp of Wettest Quarter", "Mean Temp of Driest Quarter", "Mean Temp of Warmest Quarter", "Mean Temp of Coldest Quarter", "Annual Precip", "Precip of Wettest Month", "Precip of Driest Month", "Precip Seasonality", "Precip of Wettest Quarter", "Precip of Driest Quarter", "Precip of Warmest Quarter", "Precip of Coldest Quarter") plot(myExplFuture4.5) myExplFuture4.5<-stack((myExplFuture4.5$Temp.Annual.Range/10), (myExplFuture4.5$Mean.Diurnal.Range/10), (myExplFuture4.5$Mean.Temp.of.Warmest.Quarter/10), (myExplFuture4.5$Annual.Precip/100), (myExplFuture4.5$Precip.of.Warmest.Quarter/100), dist.layer) myExplFuture4.5<-mask(myExplFuture4.5, s) myExplFuture4.5<-stack(myExplFuture4.5) names(myExplFuture4.5)<-c("ATR", "MDR", "MTWQ", "AP", "PWQ", "DIST") plot(myExplFuture4.5) myBiomodProjFuture4.5 <- BIOMOD_Projection( modeling.output = myBiomodModelOut, new.env = myExplFuture4.5, proj.name = '2050_4.5', selected.models = 'all', binary.meth = 'TSS', filtered.meth = 'TSS', compress = 'xz', clamping.mask = T, output.format = '.grd') ### Future projection plot # make some plots, sub-selected by str.grep argument plot(myBiomodProjFuture4.5, str.grep = 'GLM') plot(myBiomodProjFuture4.5, str.grep = 'GBM') plot(myBiomodProjFuture4.5, str.grep = 'RF') #2050 RCP 8.5# myExplFuture8.5.1<- getData('CMIP5', var='bio', res=5, rcp=85, model='GS', year=50) myExplFuture8.5.2<-getData('CMIP5', var='bio', res=5, rcp=85, model='HD', year=50) #myExplFuture8.5.3<- getData('CMIP5', var='bio', res=5, # rcp=85, model='GF', year=50) myExplFuture8.5.3<-stack("gf85bi501.tif", "gf85bi502.tif", "gf85bi503.tif", "gf85bi504.tif", "gf85bi505.tif", "gf85bi506.tif", "gf85bi507.tif", "gf85bi508.tif", "gf85bi509.tif", "gf85bi5010.tif", "gf85bi5011.tif", "gf85bi5012.tif", "gf85bi5013.tif", "gf85bi5014.tif", "gf85bi5015.tif", "gf85bi5016.tif", "gf85bi5017.tif", "gf85bi5018.tif", "gf85bi5019.tif") myExplFuture8.5.4<- getData('CMIP5', var='bio', res=5, rcp=85, model='IP', year=50) myExplFuture8.5.5<- getData('CMIP5', var='bio', res=5, rcp=85, model='MP', year=50) myExplFuture8.5.6<- getData('CMIP5', var='bio', res=5, rcp=85, model='NO', year=50) myExplFuture8.5.7<- getData('CMIP5', var='bio', res=5, rcp=85, model='CC', year=50) myExplFuture8.5.8<- getData('CMIP5', var='bio', res=5, rcp=85, model='BC', year=50) myExplFuture8.5<-stack(myExplFuture8.5.1, myExplFuture8.5.2, myExplFuture8.5.3, myExplFuture8.5.4, myExplFuture8.5.5, myExplFuture8.5.6, myExplFuture8.5.7, myExplFuture8.5.8) myExplFuture8.5<-subset(myExplFuture8.5, c(1,20,39,58,77,96,115,134,2,21,40,59,78,97,116,135,3,22,41,60,79,98,117,136, 4,23,42,61,80,99,118,137,5,24,43,62,81,100,119,138,6,25,44,63,82,101,120,139,7,26,45,64,83,102,121,140,8,27,46,65,84,103,122,141, 9,28,47,66,85,104,123,142,10,29,48,67,86,105,124,143,11,30,49,68,87,106,125,144,12,31,50,69,88,107,126,145,13,32,51,70,89,108,127,146, 14,33,52,71,90,109,128,147,15,34,53,72,91,110,129,148,16,35,54,73,92,111,130,149,17,36,55,74,93,112,131,150,18,37,56,75,94,113,132,151, 19,38,57,76,95,114,133,152)) myExplFuture8.5<-stack(crop(myExplFuture8.5,(extent(-180, -120, 33, 69)))) myExplFuture8.5<-stackApply(myExplFuture8.5, c(1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5, 6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,8,8,8,8,8,8,8,8,9,9,9,9,9,9,9,9,10,10,10,10,10,10,10,10, 11,11,11,11,11,11,11,11,12,12,12,12,12,12,12,12,13,13,13,13,13,13,13,13,14,14,14,14,14,14,14,14, 15,15,15,15,15,15,15,15,16,16,16,16,16,16,16,16,17,17,17,17,17,17,17,17, 18,18,18,18,18,18,18,18,19,19,19,19,19,19,19,19), fun= mean) names(myExplFuture8.5)<-c("Annual Mean Temp", "Mean Diurnal Range", "Isothermality", "Temp Seasonality", "Max Temp of Warmest", "Min Temp of Coldest", "Temp Annual Range", "Mean Temp of Wettest Quarter", "Mean Temp of Driest Quarter", "Mean Temp of Warmest Quarter", "Mean Temp of Coldest Quarter", "Annual Precip", "Precip of Wettest Month", "Precip of Driest Month", "Precip Seasonality", "Precip of Wettest Quarter", "Precip of Driest Quarter", "Precip of Warmest Quarter", "Precip of Coldest Quarter") plot(myExplFuture8.5) myExplFuture8.5<-stack((myExplFuture8.5$Temp.Annual.Range/10), (myExplFuture8.5$Mean.Diurnal.Range/10), (myExplFuture8.5$Mean.Temp.of.Warmest.Quarter/10), (myExplFuture8.5$Annual.Precip/100), (myExplFuture8.5$Precip.of.Warmest.Quarter/100), dist.layer) myExplFuture8.5<-mask(myExplFuture8.5, s) myExplFuture8.5<-stack(myExplFuture8.5) names(myExplFuture8.5)<-c("ATR", "MDR", "MTWQ", "AP", "PWQ", "DIST") plot(myExplFuture8.5) myBiomodProjFuture8.5 <- BIOMOD_Projection( modeling.output = myBiomodModelOut, new.env = myExplFuture8.5, proj.name = '2050_8.5', selected.models = 'all', binary.meth = 'TSS', filtered.meth = 'TSS', compress = 'xz', clamping.mask = T, output.format = '.grd') ### Future projection plot # make some plots, sub-selected by str.grep argument plot(myBiomodProjFuture8.5, str.grep = 'GLM') plot(myBiomodProjFuture8.5, str.grep = 'GBM') plot(myBiomodProjFuture8.5, str.grep = 'RF') ### Current period Ensemble Forecasting myBiomodEF <- BIOMOD_EnsembleForecasting( EM.output = myBiomodEM, projection.output = myBiomodProj, proj.name= 'present', binary.meth= 'TSS') ### ensemble model myBiomodEF list.files(myBiomodEF) ### plotting ensemble model # reduce layer names for plotting convegences plot(myBiomodEF) ### Future ensemble forecasting # 2050 RCP 4.5 myBiomodEF.Future4.5 <- BIOMOD_EnsembleForecasting( EM.output = myBiomodEM, projection.output = myBiomodProjFuture4.5, proj.name= '2050_4.5', binary.meth= 'TSS') plot(myBiomodEF.Future4.5) # 2050 RCP 8.5 myBiomodEF.Future8.5 <- BIOMOD_EnsembleForecasting( EM.output = myBiomodEM, projection.output = myBiomodProjFuture8.5, proj.name= '2050_8.5', binary.meth= 'TSS') plot(myBiomodEF.Future8.5) ### calculating range size changes #ensemble model currentPred<-stack('Tufted.Puffin/proj_present/proj_present_Tufted.Puffin_ensemble_TSSbin.grd') futurePred4.5<-stack('Tufted.Puffin/proj_2050_4.5/proj_2050_4.5_Tufted.Puffin_ensemble_TSSbin.grd') futurePred8.5<-stack('Tufted.Puffin/proj_2050_8.5/proj_2050_8.5_Tufted.Puffin_ensemble_TSSbin.grd') #all ensemble members individually currentPred.all<-stack("Tufted.Puffin/proj_present/proj_present_Tufted.Puffin_TSSbin.grd") futurePred.all4.5<-stack("Tufted.Puffin/proj_2050_4.5/proj_2050_4.5_Tufted.Puffin_TSSbin.grd") futurePred.all8.5<-stack("Tufted.Puffin/proj_2050_8.5/proj_2050_8.5_Tufted.Puffin_TSSbin.grd") plot(currentPred) plot(currentPred.all) plot(futurePred4.5) plot(futurePred8.5) #range change with all members included myBiomodRangeSize.all4.5<-BIOMOD_RangeSize(currentPred.all, futurePred.all4.5) myBiomodRangeSize.all8.5<-BIOMOD_RangeSize(currentPred.all, futurePred.all8.5) myBiomodRangeSize.all4.5$Compt.By.Models myBiomodRangeSize.all8.5$Compt.By.Models #range change with just ensemble model myBiomodRangeSize4.5<-BIOMOD_RangeSize(currentPred, futurePred4.5) myBiomodRangeSize8.5<-BIOMOD_RangeSize(currentPred, futurePred8.5) myBiomodRangeSize4.5$Compt.By.Models myBiomodRangeSize8.5$Compt.By.Models plot(myBiomodRangeSize4.5$Diff.By.Pixel) plot(myBiomodRangeSize8.5$Diff.By.Pixel) #California current sub section #ensemble CAC.extent.current<-stack(crop(currentPred,(extent(-126,-120,32,48.5)))) CAC.extent.future4.5<-stack(crop(futurePred4.5,(extent(-126,-120,32,48.5)))) CAC.extent.future8.5<-stack(crop(futurePred8.5,(extent(-126,-120,32,48.5)))) #individual members CAC.extent.current.all<-stack(crop(currentPred.all,(extent(-126,-120,32,48.5)))) CAC.extent.future.all4.5<-stack(crop(futurePred.all4.5,(extent(-126,-120,32,48.5)))) CAC.extent.future.all8.5<-stack(crop(futurePred.all8.5,(extent(-126,-120,32,48.5)))) Rangesize.CA4.5<-BIOMOD_RangeSize(CAC.extent.current, CAC.extent.future4.5) Rangesize.CA8.5<-BIOMOD_RangeSize(CAC.extent.current, CAC.extent.future8.5) Rangesize.CA4.5 Rangesize.CA8.5 Rangesize.CA.all4.5<-BIOMOD_RangeSize(CAC.extent.current.all, CAC.extent.future.all4.5) Rangesize.CA.all8.5<-BIOMOD_RangeSize(CAC.extent.current.all, CAC.extent.future.all8.5) #California Current sub section individual results by.model4.5<-as.data.frame(Rangesize.CA.all4.5$Compt.By.Models) by.model4.5 by.model8.5<-as.data.frame(Rangesize.CA.all8.5$Compt.By.Models) by.model8.5 #Range wide individual results by.model.all4.5<-as.data.frame(myBiomodRangeSize.all4.5$Compt.By.Models) by.model.all8.5<-as.data.frame(myBiomodRangeSize.all8.5$Compt.By.Models) #California Current sub section histogram of percent range loss #RCP 4.5 range.change4.5<-by.model4.5$PercLoss range.change4.5<-range.change4.5*-1 range.change.GLM4.5<-range.change4.5[c(1,4,7,10,13,16,19,22,25,28,31,34,37,40,43,46,49,52,55,58, 61,64,67,70,73,76,79,82,85,88,91,94,97,100,103,106,109,112, 115,118)] range.change.GBM4.5<-range.change4.5[c(2,5,8,11,14,17,20,23,26,29,32,35,38,41,44,47,50,53,56,59, 62,65,68,71,74,77,80,83,86,89,92,95,98,101,104,107,110,113, 116,119)] range.change.RF4.5<-range.change4.5[c(3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51,54,57,60, 63,66,69,72,75,78,81,84,87,90,93,96,99,102,105,108,111,114, 117,120)] range.change.total4.5<-melt(data.frame(range.change.GLM4.5, range.change.GBM4.5, range.change.RF4.5)) ggplot(range.change.total4.5, aes(x=value, fill=variable)) + geom_histogram(bins=10 )+ theme(text = element_text(size=20)) counts4.5<-table(range.change.total4.5$variable, range.change.total4.5$value) barplot(counts4.5, xlab= "% Habitat Loss", ylab= "Frequency", #main= "California Current Habitat Loss with No Migration", cex.lab=1.5, cex.axes=1.4, cex.main=1.5) ?barplot Rangesize.CA4.5$Compt.By.Models Rangesize.CA.all4.5$Compt.By.Models #RCP 8.5 range.change8.5<-by.model8.5$PercLoss range.change8.5<-range.change8.5*-1 range.change.GLM8.5<-range.change8.5[c(1,4,7,10,13,16,19,22,25,28,31,34,37,40,43,46,49,52,55,58, 61,64,67,70,73,76,79,82,85,88,91,94,97,100,103,106,109,112, 115,118)] range.change.GBM8.5<-range.change8.5[c(2,5,8,11,14,17,20,23,26,29,32,35,38,41,44,47,50,53,56,59, 62,65,68,71,74,77,80,83,86,89,92,95,98,101,104,107,110,113, 116,119)] range.change.RF8.5<-range.change8.5[c(3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51,54,57,60, 63,66,69,72,75,78,81,84,87,90,93,96,99,102,105,108,111,114, 117,120)] range.change.total8.5<-melt(data.frame(range.change.GLM8.5, range.change.GBM8.5, range.change.RF8.5)) ggplot(range.change.total8.5, aes(x=value, fill=variable)) + geom_histogram(bins=8, colour= "black")+ scale_fill_discrete(name="Model\nAlgorithm", labels=c("GLM", "GBM", "RF"))+ labs(x="% Habitat Loss", y="Frequency")+theme(text = element_text(size=15)) counts8.5<-table(range.change.total8.5$variable, range.change.total8.5$value) hist(range.change.total8.5$value, xlab= "% Habitat Loss", ylab= "Frequency", main= "California Current Habitat Loss with No Migration", breaks=10, col=range.change.total8.5$variable) Rangesize.CA8.5$Compt.By.Models Rangesize.CA.all8.5$Compt.By.Models #Plot of projected loss on California Current #RCP 4.5 plot(Rangesize.CA4.5$Diff.By.Pixel) plot(Rangesize.CA.all4.5$Diff.By.Pixel$layer.23) plot(Rangesize.CA4.5$Diff.By.Pixel, col=c('black', 'magenta', 'grey')) sp::spplot(Rangesize.CA4.5$Diff.By.Pixel, key.space= "right", scales=list(draw=T), col.regions= c("Black", "Green", "Grey"), cuts=2, main= "RCP 4.5") #RCP 8.5 plot(Rangesize.CA8.5$Diff.By.Pixel) plot(Rangesize.CA.all8.5$Diff.By.Pixel$layer.23) plot(Rangesize.CA8.5$Diff.By.Pixel, col=c('black', 'magenta', 'grey')) sp::spplot(Rangesize.CA8.5$Diff.By.Pixel, key.space= "right", scales=list(draw=T), col.regions= c("Black", "Green", "Grey"), cuts=2) # % habitat loss #for individual ensemble members only #California Current sub section analysis #RCP 4.5 GLM.loss.avg4.5<-mean(by.model4.5[c(1,4,7,10,13,16,19,22,25,28,31,34,37,40,43,46,49,52,55,58, 61,64,67,70,73,76,79,82,85,88,91,94,97,100,103,106,109,112, 115,118),7]) GLM.loss.std4.5<-sd(by.model4.5[c(1,4,7,10,13,16,19,22,25,28,31,34,37,40,43,46,49,52,55,58, 61,64,67,70,73,76,79,82,85,88,91,94,97,100,103,106,109,112, 115,118),7]) GBM.loss.avg4.5<-(by.model4.5[c(2,5,8,11,14,17,20,23,26,29,32,35,38,41,44,47,50,53,56,59, 62,65,68,71,74,77,80,83,86,89,92,95,98,101,104,107,110,113, 116),7]) GBM.loss.avg4.5<-mean(GBM.loss.avg4.5, na.rm=T) GBM.loss.std4.5<-(by.model4.5[c(2,5,8,11,14,17,20,23,26,29,32,35,38,41,44,47,50,53,56,59, 62,65,68,71,74,77,80,83,86,89,92,95,98,101,104,107,110,113, 116),7]) GBM.loss.std4.5<-sd(GBM.loss.std4.5, na.rm=T) RF.loss.avg4.5<-mean(by.model4.5[c(3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51,54,57,60, 63,66,69,72,75,78,81,84,87,90,93,96,99,102,105,108,111,114, 117,120),7]) RF.loss.std4.5<-sd(by.model4.5[c(3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51,54,57,60, 63,66,69,72,75,78,81,84,87,90,93,96,99,102,105,108,111,114, 117,120),7]) #RCP 8.5 GLM.loss.avg8.5<-mean(by.model8.5[c(1,4,7,10,13,16,19,22,25,28,31,34,37,40,43,46,49,52,55,58, 61,64,67,70,73,76,79,82,85,88,91,94,97,100,103,106,109,112, 115,118),7]) GLM.loss.std8.5<-sd(by.model8.5[c(1,4,7,10,13,16,19,22,25,28,31,34,37,40,43,46,49,52,55,58, 61,64,67,70,73,76,79,82,85,88,91,94,97,100,103,106,109,112, 115,118),7]) GBM.loss.avg8.5<-(by.model8.5[c(2,5,8,11,14,17,20,23,26,29,32,35,38,41,44,47,50,53,56,59, 62,65,68,71,74,77,80,83,86,89,92,95,98,101,104,107,110,113, 116),7]) GBM.loss.avg8.5<-mean(GBM.loss.avg8.5, na.rm=T) GBM.loss.std8.5<-(by.model8.5[c(2,5,8,11,14,17,20,23,26,29,32,35,38,41,44,47,50,53,56,59, 62,65,68,71,74,77,80,83,86,89,92,95,98,101,104,107,110,113, 116),7]) GBM.loss.std8.5<-sd(GBM.loss.std8.5, na.rm=T) RF.loss.avg8.5<-mean(by.model8.5[c(3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51,54,57,60, 63,66,69,72,75,78,81,84,87,90,93,96,99,102,105,108,111,114, 117,120),7]) RF.loss.std8.5<-sd(by.model8.5[c(3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51,54,57,60, 63,66,69,72,75,78,81,84,87,90,93,96,99,102,105,108,111,114, 117,120),7]) #Range wide analysis # RCP 4.5 GLM.all.loss.avg4.5<-mean(by.model.all4.5[c(1,4,7,10,13,16,19,22,25,28,31,34,37,40,43,46,49,52,55,58, 61,64,67,70,73,76,79,82,85,88,91,94,97,100,103,106,109,112, 115,118),7]) GLM.all.loss.std4.5<-sd(by.model.all4.5[c(1,4,7,10,13,16,19,22,25,28,31,34,37,40,43,46,49,52,55,58, 61,64,67,70,73,76,79,82,85,88,91,94,97,100,103,106,109,112, 115,118),7]) GBM.all.loss.avg4.5<-mean(by.model.all4.5[c(2,5,8,11,14,17,20,23,26,29,32,35,38,41,44,47,50,53,56,59, 62,65,68,71,74,77,80,83,86,89,92,95,98,101,104,107,110,113, 116,119),7]) GBM.all.loss.std4.5<-sd(by.model.all4.5[c(2,5,8,11,14,17,20,23,26,29,32,35,38,41,44,47,50,53,56,59, 62,65,68,71,74,77,80,83,86,89,92,95,98,101,104,107,110,113, 116,119),7]) RF.all.loss.avg4.5<-mean(by.model.all4.5[c(3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51,54,57,60, 63,66,69,72,75,78,81,84,87,90,93,96,99,102,105,108,111,114, 117,120),7]) RF.all.loss.std4.5<-sd(by.model.all4.5[c(3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51,54,57,60, 63,66,69,72,75,78,81,84,87,90,93,96,99,102,105,108,111,114, 117,120),7]) # RCP 8.5 GLM.all.loss.avg8.5<-mean(by.model.all8.5[c(1,4,7,10,13,16,19,22,25,28,31,34,37,40,43,46,49,52,55,58, 61,64,67,70,73,76,79,82,85,88,91,94,97,100,103,106,109,112, 115,118),7]) GLM.all.loss.std8.5<-sd(by.model.all8.5[c(1,4,7,10,13,16,19,22,25,28,31,34,37,40,43,46,49,52,55,58, 61,64,67,70,73,76,79,82,85,88,91,94,97,100,103,106,109,112, 115,118),7]) GBM.all.loss.avg8.5<-mean(by.model.all8.5[c(2,5,8,11,14,17,20,23,26,29,32,35,38,41,44,47,50,53,56,59, 62,65,68,71,74,77,80,83,86,89,92,95,98,101,104,107,110,113, 116,119),7]) GBM.all.loss.std8.5<-sd(by.model.all8.5[c(2,5,8,11,14,17,20,23,26,29,32,35,38,41,44,47,50,53,56,59, 62,65,68,71,74,77,80,83,86,89,92,95,98,101,104,107,110,113, 116,119),7]) RF.all.loss.avg8.5<-mean(by.model.all8.5[c(3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51,54,57,60, 63,66,69,72,75,78,81,84,87,90,93,96,99,102,105,108,111,114, 117,120),7]) RF.all.loss.std8.5<-sd(by.model.all8.5[c(3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51,54,57,60, 63,66,69,72,75,78,81,84,87,90,93,96,99,102,105,108,111,114, 117,120),7]) #Google Maps library(raster); library(ggplot2); library(ggmap) #devtools::install_github("dkahle/ggmap") # CC_Present CC_present<-raster("Tufted.Puffin/proj_present/proj_present_Tufted.Puffin_ensemble.gri") plot(CC_present) CC_present CC.extent<-extent(-125, -120, 34, 50) #create cropping extent CC_present<-crop(CC_present, CC.extent) #crop to CA current plot(CC_present) #create data frame from raster data; note we could rewrite this easily to pull the parameters from the crop extent (y), if we want to highlight different spatial areas... df.CC.present<-data.frame(lat=rep(seq(34,50, length.out=192), each=60), long=rep(seq(-120, -125, length.out=60), times=192), value=rev(as.vector(CC_present)) ) df.CC.present<-df.CC.present[!is.na(df.CC.present$value),] #remove missing values df.CC.present<-df.CC.present[df.CC.present$value>100,] #remove very low values df.CC.present$value<-cut(df.CC.present$value, 4) #bin data into 4 bins reflecting likelihood of occupancy #create color palette plotcolors <- colorRampPalette(c("indianred3", "orange", "yellow", "green"))(4) #now, to put this on a google map map.CC <- get_map(location = c(lon = mean(df.CC.present$long), lat = mean(df.CC.present$lat)), zoom = 5, maptype = "terrain-background", source = "stamen" ) plotmap.CC.present<-ggmap(map.CC, extent="panel") plotmap.CC.present<-plotmap.CC.present+ geom_point(data=df.CC.present, aes(x=long, y=lat, color=value))+ annotate("text",x=-132, y=50,label="A",color="Black", size=15)+ annotate("text", x=-118.5, y=33.25, label="© Stamen Design, under CC BY 3.0. Data © OpenStreetMap, under ODbL.", size=2.0)+ scale_color_manual(values= plotcolors, name="Probability \nof Occupancy", labels=c("10 - 27%","27% - 44%","44 - 61%","61 - 78%")) plotmap.CC.present ggsave("CalCurrent_present_plot_ggmaps.pdf") # CC-2050-4.5 CC_2050_4.5<-raster("Tufted.Puffin/proj_2050_4.5/proj_2050_4.5_Tufted.Puffin_ensemble.gri") plot(CC_2050_4.5) CC_2050_4.5 CC_2050_4.5<-crop(CC_2050_4.5, CC.extent) #crop to CA current plot(CC_2050_4.5) #create data frame from raster data df.CC.2050.4.5<-data.frame(lat=rep(seq(34,50, length.out=192), each=60), long=rep(seq(-120, -125, length.out=60), times=192), value=rev(as.vector(CC_2050_4.5)) ) df.CC.2050.4.5<-df.CC.2050.4.5[!is.na(df.CC.2050.4.5$value),] #remove missing values df.CC.2050.4.5<-df.CC.2050.4.5[df.CC.2050.4.5$value>100,] #remove very low values df.CC.2050.4.5$value<-cut(df.CC.2050.4.5$value, c(100,272,442,612)) #df.CC.2050.4.5$value<-cut(df.CC.2050.4.5$value, 4) #bin data into 4 bins reflecting likelihood of occupancy #create color palette plotcolors <- colorRampPalette(c("indianred3", "orange", "yellow", "green"))(4) #now, to put this on a google map map.CC.2050.4.5 <- get_map(location = c(lon = mean(df.CC.2050.4.5$long), lat = mean(df.CC.2050.4.5$lat)), zoom = 5, maptype = "terrain-background", source = "stamen" ) plotmap.CC.2050.4.5<-ggmap(map.CC.2050.4.5) plotmap.CC.2050.4.5<-plotmap.CC.2050.4.5+ geom_point(data=df.CC.2050.4.5, aes(x=long, y=lat, color=value))+ annotate("text",x=-132, y=50,label="B",color="Black", size=15)+ annotate("text", x=-118.5, y=33.25, label="© Stamen Design, under CC BY 3.0. Data © OpenStreetMap, under ODbL.", size=2.0)+ scale_color_manual(values= plotcolors, name="Probability \nof Occupancy", labels=c("10 - 27%","27% - 44%","44 - 61%","61 - 78%")) plotmap.CC.2050.4.5 ggsave("CalCurrent_2050.4.5_plot_ggmaps.pdf") # CC-2050-8.5 CC_2050_8.5<-raster("Tufted.Puffin/proj_2050_8.5/proj_2050_8.5_Tufted.Puffin_ensemble.gri") plot(CC_2050_8.5) CC_2050_8.5 CC_2050_8.5<-crop(CC_2050_8.5, CC.extent) #crop to CA current plot(CC_2050_8.5) #create data frame from raster data df.CC.2050.8.5<-data.frame(lat=rep(seq(34,50, length.out=192), each=60), long=rep(seq(-120, -125, length.out=60), times=192), value=rev(as.vector(CC_2050_8.5)) ) df.CC.2050.8.5<-df.CC.2050.8.5[!is.na(df.CC.2050.8.5$value),] #remove missing values df.CC.2050.8.5<-df.CC.2050.8.5[df.CC.2050.8.5$value>100,] #remove very low values df.CC.2050.8.5$value<-cut(df.CC.2050.8.5$value, c(100,272,442,612)) #bin data into 4 bins reflecting likelihood of occupancy #create color palette plotcolors <- colorRampPalette(c("indianred3", "orange", "yellow", "green"))(4) #now, to put this on a google map map.CC.2050.8.5 <- get_map(location = c(lon = mean(df.CC.2050.8.5$long), lat = mean(df.CC.2050.8.5$lat)), zoom = 5, maptype = "terrain-background", source = "stamen" ) plotmap.CC.2050.8.5<-ggmap(map.CC.2050.8.5) plotmap.CC.2050.8.5<-plotmap.CC.2050.8.5+ geom_point(data=df.CC.2050.8.5, aes(x=long, y=lat, color=value))+ annotate("text",x=-132, y=50,label="C",color="Black", size=15)+ annotate("text", x=-118.5, y=33.25, label="© Stamen Design, under CC BY 3.0. Data © OpenStreetMap, under ODbL.", size=2.0)+ scale_color_manual(values= plotcolors, name="Probability \nof Occupancy", labels=c("10 - 27%","27% - 44%","44 - 61%","61 - 78%")) plotmap.CC.2050.8.5 ggsave("CalCurrent_2050.8.5_plot_ggmaps.pdf") # Multiple plot function # # ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects) # - cols: Number of columns in layout # - layout: A matrix specifying the layout. If present, 'cols' is ignored. # # If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE), # then plot 1 will go in the upper left, 2 will go in the upper right, and # 3 will go all the way across the bottom. # multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) { library(grid) # Make a list from the ... arguments and plotlist plots <- c(list(...), plotlist) numPlots = length(plots) # If layout is NULL, then use 'cols' to determine layout if (is.null(layout)) { # Make the panel # ncol: Number of columns of plots # nrow: Number of rows needed, calculated from # of cols layout <- matrix(seq(1, cols * ceiling(numPlots/cols)), ncol = cols, nrow = ceiling(numPlots/cols)) } if (numPlots==1) { print(plots[[1]]) } else { # Set up the page grid.newpage() pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout)))) # Make each plot, in the correct location for (i in 1:numPlots) { # Get the i,j matrix positions of the regions that contain this subplot matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE)) print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row, layout.pos.col = matchidx$col)) } } } multiplot(plotmap.CC.present, plotmap.CC.2050.4.5, plotmap.CC.2050.8.5, cols=2, rows=2) #Range Wide Maps ## Range Wide-Present RW.current<-raster("Tufted.Puffin/proj_present/proj_present_Tufted.Puffin_ensemble.gri") RW.current plot(RW.current) RW.extent<-extent(-180, -120, 33, 68) #create cropping extent RW.current<-crop(RW.current, RW.extent) #crop to RW plot(RW.current) #create data frame from raster data; note we could rewrite this easily to pull the parameters from the crop extent (y), if we want to highlight different spatial areas... Range.data<-data.frame(lat=rep(seq(33,68, length.out=420), each=720), long=rep(seq(-120, -180, length.out=720), times=420), value=rev(as.vector(RW.current)) ) Range.data<-Range.data[complete.cases(Range.data),] Range.data<-Range.data[Range.data$value>100,] #remove very low values Range.data$value<-cut(Range.data$value, 4) #bin data into 4 bins reflecting likelihood of occupancy #create color palette plotcolors <- colorRampPalette(c("orange", "yellow", "green", "dark green"))(4) #now, to put this on a google map RW.map<-get_googlemap(center= c(lon=mean(Range.data$long), lat=mean(Range.data$lat-3)), zoom=3, size=c(500, 400), scale=2, maptype = "terrain") plotmap.RW.present<-ggmap(RW.map, extent="panel") plotmap.RW.present<-plotmap.RW.present+ geom_point(data=Range.data, aes(x=long, y=lat, color=value))+ annotate("text",x=-187, y=66,label="A",color="Black", size=15)+ scale_color_manual(values= plotcolors, name="Probability \nof Occupancy", labels=c("10 - 32%","32% - 54%","54 - 76%","76 - 98%")) plotmap.RW.present ggsave("Rangewide_present_plot_ggmaps.pdf") ## Range Wide-4.5 RW.4.5<-raster("Tufted.Puffin/proj_2050_4.5/proj_2050_4.5_Tufted.Puffin_ensemble.gri") RW.4.5 plot(RW.4.5) RW.4.5<-crop(RW.4.5, RW.extent) #crop to RW plot(RW.4.5) #create data frame from raster data Range.data.4.5<-data.frame(lat=rep(seq(33,68, length.out=420), each=720), long=rep(seq(-120, -180, length.out=720), times=420), value=rev(as.vector(RW.4.5)) ) Range.data.4.5<-Range.data.4.5[complete.cases(Range.data.4.5),] Range.data.4.5<-Range.data.4.5[Range.data.4.5$value>100,] #remove very low values Range.data.4.5$value<-cut(Range.data.4.5$value, c(100,320,540,759,979)) #bin data into 4 bins reflecting likelihood of occupancy #create color palette plotcolors <- colorRampPalette(c("orange", "yellow", "green", "dark green"))(4) #now, to put this on a google map RW.map.4.5<-get_googlemap(center= c(lon=mean(Range.data.4.5$long), lat=mean(Range.data.4.5$lat-3)), zoom=3, size=c(500, 400), scale=2, maptype = "terrain") plotmap.RW.4.5<-ggmap(RW.map.4.5, extent="panel") plotmap.RW.4.5<-plotmap.RW.4.5+ geom_point(data=Range.data.4.5, aes(x=long, y=lat, color=value))+ annotate("text",x=-190, y=66,label="B",color="Black", size=15)+ scale_color_manual(values= plotcolors, name="Probability \nof Occupancy", labels=c("10 - 32%","32% - 54%","54 - 76%","76 - 98%")) plotmap.RW.4.5 ggsave("Rangewide_4.5_plot_ggmaps.pdf") ## Range Wide-8.5 RW.8.5<-raster("Tufted.Puffin/proj_2050_8.5/proj_2050_8.5_Tufted.Puffin_ensemble.gri") RW.8.5 plot(RW.8.5) RW.8.5<-crop(RW.8.5, RW.extent) #crop to RW plot(RW.8.5) #create data frame from raster data Range.data.8.5<-data.frame(lat=rep(seq(33,68, length.out=420), each=720), long=rep(seq(-120, -180, length.out=720), times=420), value=rev(as.vector(RW.8.5)) ) Range.data.8.5<-Range.data.8.5[complete.cases(Range.data.8.5),] Range.data.8.5<-Range.data.8.5[Range.data.8.5$value>100,] #remove very low values Range.data.8.5$value<-cut(Range.data.8.5$value, c(100,320,540,759,979)) #bin data into 4 bins reflecting likelihood of occupancy #create color palette plotcolors <- colorRampPalette(c("orange", "yellow", "green", "dark green"))(4) #now, to put this on a google map RW.map.8.5<-get_googlemap(center= c(lon=mean(Range.data.8.5$long), lat=mean(Range.data.8.5$lat-3)), zoom=3, size=c(500, 400), scale=2, maptype = "terrain") plotmap.RW.8.5<-ggmap(RW.map.8.5, extent="panel") plotmap.RW.8.5<-plotmap.RW.8.5+ geom_point(data=Range.data.8.5, aes(x=long, y=lat, color=value))+ annotate("text",x=-190, y=66,label="C",color="Black", size=15)+ scale_color_manual(values= plotcolors, name="Probability \nof Occupancy", labels=c("10 - 32%","32% - 54%","54 - 76%","76 - 98%")) plotmap.RW.8.5 ggsave("Rangewide_8.5_plot_ggmaps.pdf") multiplot(plotmap.RW.present, plotmap.RW.4.5, plotmap.RW.8.5, cols=2, rows=2) #library(gridExtra) #png(file="mrdq.png",w=1800,h=1800, res=300) #grid.newpage() #v1<-viewport(width = 1, height = 1, x = 0.5, y = 0.5) #plot area for the main map #v2<-viewport(width = 0.35, height = 0.35, x = 0.1, y = 0.1) #plot area for the inset map #print(plotmap.CC.present,vp=v1) #print(plotmap.RW.present,vp=v2) #dev.off() # PAST CLIMATE RUN, Hindcast #Mean temp of warmest quarter pastExpl.twarm1<-stack("cru_ts4.01.1911.1920.tmp.dat.nc", varname="tmp") pastExpl.twarm1<-stack(crop(pastExpl.twarm1,(extent(-180,-120,32,69)))) pastExpl.twarm2<-stack("cru_ts4.01.1921.1930.tmp.dat.nc",varname="tmp") pastExpl.twarm2<-stack(crop(pastExpl.twarm2,(extent(-180,-120,32,69)))) pastExpl.twarm3<-stack("cru_ts4.01.1931.1940.tmp.dat.nc",varname="tmp") pastExpl.twarm3<-stack(crop(pastExpl.twarm3,(extent(-180,-120,32,69)))) pastExpl.twarm4<-stack("cru_ts4.01.1941.1950.tmp.dat.nc",varname="tmp") pastExpl.twarm4<-stack(crop(pastExpl.twarm4,(extent(-180,-120,32,69)))) indices.past.twarm <- format(as.Date(names(pastExpl.twarm1), format = "X%Y.%m.%d"), format = "%m") indices.past.twarm <- as.numeric(indices.past.twarm) pastExpl.twarm1<- stackApply(pastExpl.twarm1, indices.past.twarm, fun = mean) pastExpl.twarm1<-stack(pastExpl.twarm1$index_6, pastExpl.twarm1$index_7, pastExpl.twarm1$index_8) pastExpl.twarm2<- stackApply(pastExpl.twarm2, indices.past.twarm, fun = mean) pastExpl.twarm2<-stack(pastExpl.twarm2$index_6, pastExpl.twarm2$index_7, pastExpl.twarm2$index_8) pastExpl.twarm3<- stackApply(pastExpl.twarm3, indices.past.twarm, fun = mean) pastExpl.twarm3<-stack(pastExpl.twarm3$index_6, pastExpl.twarm3$index_7, pastExpl.twarm3$index_8) pastExpl.twarm4<- stackApply(pastExpl.twarm4, indices.past.twarm, fun = mean) pastExpl.twarm4<-stack(pastExpl.twarm4$index_6, pastExpl.twarm4$index_7, pastExpl.twarm4$index_8) pastExpltwarm<-stack(pastExpl.twarm1, pastExpl.twarm2, pastExpl.twarm3, pastExpl.twarm4) pastExpltwarm<-stackApply(pastExpltwarm, indices = 12, fun = mean) plot(pastExpltwarm) #Precip of warmest quarter pastExpl.pwarm1<-stack("cru_ts4.01.1911.1920.pre.dat.nc",varname="pre") pastExpl.pwarm1<-stack(crop(pastExpl.pwarm1,(extent(-180,-120,32,69)))) pastExpl.pwarm2<-stack("cru_ts4.01.1921.1930.pre.dat.nc",varname="pre") pastExpl.pwarm2<-stack(crop(pastExpl.pwarm2,(extent(-180,-120,32,69)))) pastExpl.pwarm3<-stack("cru_ts4.01.1931.1940.pre.dat.nc",varname="pre") pastExpl.pwarm3<-stack(crop(pastExpl.pwarm3,(extent(-180,-120,32,69)))) pastExpl.pwarm4<-stack("cru_ts4.01.1941.1950.pre.dat.nc",varname="pre") pastExpl.pwarm4<-stack(crop(pastExpl.pwarm4,(extent(-180,-120,32,69)))) indices.past.pwarm <- format(as.Date(names(pastExpl.pwarm1), format = "X%Y.%m.%d"), format = "%m") indices.past.pwarm <- as.numeric(indices.past.pwarm) pastExpl.pwarm1<- stackApply(pastExpl.pwarm1, indices.past.pwarm, fun = mean) pastExpl.pwarm1<-stack(pastExpl.pwarm1$index_6, pastExpl.pwarm1$index_7, pastExpl.pwarm1$index_8) pastExpl.pwarm2<- stackApply(pastExpl.pwarm2, indices.past.pwarm, fun = mean) pastExpl.pwarm2<-stack(pastExpl.pwarm2$index_6, pastExpl.pwarm2$index_7, pastExpl.pwarm2$index_8) pastExpl.pwarm3<- stackApply(pastExpl.pwarm3, indices.past.pwarm, fun = mean) pastExpl.pwarm3<-stack(pastExpl.pwarm3$index_6, pastExpl.pwarm3$index_7, pastExpl.pwarm3$index_8) pastExpl.pwarm4<- stackApply(pastExpl.pwarm4, indices.past.pwarm, fun = mean) pastExpl.pwarm4<-stack(pastExpl.pwarm4$index_6, pastExpl.pwarm4$index_7, pastExpl.pwarm4$index_8) pastExplpwarm<-stack(pastExpl.pwarm1, pastExpl.pwarm2, pastExpl.pwarm3, pastExpl.pwarm4) pastExplpwarm<-stackApply(pastExplpwarm, indices = 12, fun = mean) plot(pastExplpwarm) #annual precip pastExpl.annp1<-stack("cru_ts4.01.1911.1920.pre.dat.nc",varname="pre") pastExpl.annp1<-stack(crop(pastExpl.annp1,(extent(-180,-120,32,69)))) pastExpl.annp2<-stack("cru_ts4.01.1921.1930.pre.dat.nc",varname="pre") pastExpl.annp2<-stack(crop(pastExpl.annp2,(extent(-180,-120,32,69)))) pastExpl.annp3<-stack("cru_ts4.01.1931.1940.pre.dat.nc",varname="pre") pastExpl.annp3<-stack(crop(pastExpl.annp3,(extent(-180,-120,32,69)))) pastExpl.annp4<-stack("cru_ts4.01.1941.1950.pre.dat.nc",varname="pre") pastExpl.annp4<-stack(crop(pastExpl.annp4,(extent(-180,-120,32,69)))) indices.past.annp <- format(as.Date(names(pastExpl.annp1), format = "X%Y.%m.%d"), format = "%m") indices.past.annp <- as.numeric(indices.past.annp) pastExpl.annp1<- stackApply(pastExpl.annp1, indices.past.pwarm, fun = mean) pastExpl.annp2<- stackApply(pastExpl.annp2, indices.past.pwarm, fun = mean) pastExpl.annp3<- stackApply(pastExpl.annp3, indices.past.pwarm, fun = mean) pastExpl.annp4<- stackApply(pastExpl.annp4, indices.past.pwarm, fun = mean) pastExplannp<-stack(pastExpl.annp1, pastExpl.annp2, pastExpl.annp3, pastExpl.annp4) pastExplannp<-stackApply(pastExplannp, indices = 12, fun = mean) plot(pastExplannp) #mean diurnal range pastExpl.dtr1<-stack("cru_ts4.01.1911.1920.dtr.dat.nc",varname="dtr") pastExpl.dtr1<-stack(crop(pastExpl.dtr1,(extent(-180,-120,32,69)))) pastExpl.dtr2<-stack("cru_ts4.01.1921.1930.dtr.dat.nc",varname="dtr") pastExpl.dtr2<-stack(crop(pastExpl.dtr2,(extent(-180,-120,32,69)))) pastExpl.dtr3<-stack("cru_ts4.01.1931.1940.dtr.dat.nc",varname="dtr") pastExpl.dtr3<-stack(crop(pastExpl.dtr3,(extent(-180,-120,32,69)))) pastExpl.dtr4<-stack("cru_ts4.01.1941.1950.dtr.dat.nc",varname="dtr") pastExpl.dtr4<-stack(crop(pastExpl.dtr4,(extent(-180,-120,32,69)))) indices.past.dtr <- format(as.Date(names(pastExpl.dtr1), format = "X%Y.%m.%d"), format = "%m") indices.past.dtr <- as.numeric(indices.past.dtr) pastExpl.dtr1<- stackApply(pastExpl.dtr1, indices.past.dtr, fun = mean) pastExpl.dtr2<- stackApply(pastExpl.dtr2, indices.past.dtr, fun = mean) pastExpl.dtr3<- stackApply(pastExpl.dtr3, indices.past.dtr, fun = mean) pastExpl.dtr4<- stackApply(pastExpl.dtr4, indices.past.dtr, fun = mean) pastExpldtr<-stack(pastExpl.dtr1, pastExpl.dtr2, pastExpl.dtr3, pastExpl.dtr4) pastExpldtr<-stackApply(pastExpldtr, indices = 12, fun = mean) plot(pastExpldtr) #temperature annual range pastExpl.tmax1<-stack("cru_ts4.01.1911.1920.tmx.dat.nc", varname="tmx") pastExpl.tmax1<-stack(crop(pastExpl.tmax1,(extent(-180,-120,32,69)))) pastExpl.tmax2<-stack("cru_ts4.01.1921.1930.tmx.dat.nc", varname="tmx") pastExpl.tmax2<-stack(crop(pastExpl.tmax2,(extent(-180,-120,32,69)))) pastExpl.tmax3<-stack("cru_ts4.01.1931.1940.tmx.dat.nc", varname="tmx") pastExpl.tmax3<-stack(crop(pastExpl.tmax3,(extent(-180,-120,32,69)))) pastExpl.tmax4<-stack("cru_ts4.01.1941.1950.tmx.dat.nc", varname="tmx") pastExpl.tmax4<-stack(crop(pastExpl.tmax4,(extent(-180,-120,32,69)))) indices.past.tmax <- format(as.Date(names(pastExpl.tmax1), format = "X%Y.%m.%d"), format = "X%Y") pastExpl.tmax1<- stackApply(pastExpl.tmax1, indices.past.tmax, fun = max) pastExpl.tmax2<- stackApply(pastExpl.tmax2, indices.past.tmax, fun = max) pastExpl.tmax3<- stackApply(pastExpl.tmax3, indices.past.tmax, fun = max) pastExpl.tmax4<- stackApply(pastExpl.tmax4, indices.past.tmax, fun = max) pastExpl.tmin1<-stack("cru_ts4.01.1911.1920.tmn.dat.nc", varname="tmn") pastExpl.tmin1<-stack(crop(pastExpl.tmin1,(extent(-180,-120,32,69)))) pastExpl.tmin2<-stack("cru_ts4.01.1921.1930.tmn.dat.nc", varname="tmn") pastExpl.tmin2<-stack(crop(pastExpl.tmin2,(extent(-180,-120,32,69)))) pastExpl.tmin3<-stack("cru_ts4.01.1931.1940.tmn.dat.nc", varname="tmn") pastExpl.tmin3<-stack(crop(pastExpl.tmin3,(extent(-180,-120,32,69)))) pastExpl.tmin4<-stack("cru_ts4.01.1941.1950.tmn.dat.nc", varname="tmn") pastExpl.tmin4<-stack(crop(pastExpl.tmin4,(extent(-180,-120,32,69)))) indices.past.tmin <- format(as.Date(names(pastExpl.tmin1), format = "X%Y.%m.%d"), format = "X%Y") pastExpl.tmin1<- stackApply(pastExpl.tmin1, indices.past.tmin, fun = min) pastExpl.tmin2<- stackApply(pastExpl.tmin2, indices.past.tmin, fun = min) pastExpl.tmin3<- stackApply(pastExpl.tmin3, indices.past.tmin, fun = min) pastExpl.tmin4<- stackApply(pastExpl.tmin4, indices.past.tmin, fun = min) TAR1<-overlay(pastExpl.tmax1, pastExpl.tmin1, fun=function(x,y){return(x-y)}) TAR2<-overlay(pastExpl.tmax2, pastExpl.tmin2, fun=function(x,y){return(x-y)}) TAR3<-overlay(pastExpl.tmax3, pastExpl.tmin3, fun=function(x,y){return(x-y)}) TAR4<-overlay(pastExpl.tmax4, pastExpl.tmin4, fun=function(x,y){return(x-y)}) pastExplTAR<-stack(TAR1, TAR2, TAR3, TAR4) pastExplTAR<-stackApply(pastExplTAR, indices = 10, fun = mean) plot(pastExplTAR) library(spatial.tools) pastExpl<-stack(pastExplTAR, pastExpldtr, pastExpltwarm, pastExplannp, pastExplpwarm) pastExpl<-merged.layers<-spatial_sync_raster(pastExpl, myExpl, method= "ngb") pastExpl<-stack(pastExpl, dist.layer) names(pastExpl)<-c("ATR", "MDR", "MTWQ", "AP", "PWQ", "DIST") plot(pastExpl) pastExpl<-stack(pastExpl$ATR, pastExpl$MDR, pastExpl$MTWQ, pastExpl$AP/10, pastExpl$PWQ/10, dist.layer) pastExpl<-mask(pastExpl, s) pastExpl<-stack(pastExpl) names(pastExpl)<-c("ATR", "MDR", "MTWQ", "AP", "PWQ", "DIST") plot(pastExpl) #Past model by model projection myBiomodProjPast <- BIOMOD_Projection( modeling.output = myBiomodModelOut, new.env = pastExpl, proj.name = 'past', selected.models = 'all', binary.meth = 'TSS', filtered.meth = 'TSS', compress = 'xz', clamping.mask = T, output.format = '.grd') #Past ensemble projection myBiomodEF.Past <- BIOMOD_EnsembleForecasting( EM.output = myBiomodEM, projection.output = myBiomodProjPast, proj.name= 'past', binary.meth= 'TSS') plot(myBiomodEF.Past) #ensemble model currentPred.comp<-stack('Tufted.Puffin/proj_present/proj_present_Tufted.Puffin_ensemble_TSSbin.grd') pastPred<-stack('Tufted.Puffin/proj_past/proj_past_Tufted.Puffin_ensemble_TSSbin.grd') #all ensemble members individually currentPred.all.comp<-stack("Tufted.Puffin/proj_present/proj_present_Tufted.Puffin_TSSbin.grd") pastPred.all<-stack("Tufted.Puffin/proj_past/proj_past_Tufted.Puffin_TSSbin.grd") #range change with all members included myBiomodRangeSizepast.all<-BIOMOD_RangeSize(currentPred.all.comp, pastPred.all) myBiomodRangeSizepast.all$Compt.By.Models #range change with just ensemble model myBiomodRangeSizepast<-BIOMOD_RangeSize(currentPred.comp, pastPred) myBiomodRangeSizepast$Compt.By.Models plot(myBiomodRangeSizepast$Diff.By.Pixel) #California current sub section #ensemble CAC.extent.current.comp<-stack(crop(currentPred.comp,(extent(-126,-120,33,50)))) CAC.extent.past<-stack(crop(pastPred,(extent(-126,-120,33,50)))) #individual members CAC.extent.current.all.comp<-stack(crop(currentPred.all.comp,(extent(-126,-120,33,50)))) CAC.extent.past.all<-stack(crop(pastPred.all,(extent(-126,-120,33,50)))) Rangesize.CApast<-BIOMOD_RangeSize(CAC.extent.current.comp, CAC.extent.past) Rangesize.CA.allpast<-BIOMOD_RangeSize(CAC.extent.current.all.comp, CAC.extent.past.all) sp::spplot(Rangesize.CApast$Diff.By.Pixel, key.space= "right", scales=list(draw=T), col.regions= c("Black", "Green", "Grey", "Blue"), cuts=3, main= "1950 Climate") # CC_Past CC_past<-raster("Tufted.Puffin/proj_past/proj_past_Tufted.Puffin_ensemble.gri") CC.past<-spatial_sync_raster(CC_past, CC_present) plot(CC_past) CC_past CC.extent<-extent(-125, -120, 34, 50) #create cropping extent CC_past<-crop(CC_past, CC.extent) #crop to CA current plot(CC_past) #create data frame from raster data; note we could rewrite this easily to pull the parameters from the crop extent (y), if we want to highlight different spatial areas... df.CC.past<-data.frame(lat=rep(seq(34,50, length.out=192), each=60), long=rep(seq(-120, -125, length.out=60), times=192), value=rev(as.vector(CC_past)) ) df.CC.past<-df.CC.past[!is.na(df.CC.past$value),] #remove missing values df.CC.past<-df.CC.past[df.CC.past$value>100,] #remove very low values df.CC.past$value<-cut(df.CC.past$value, c(100,272,442,612,782)) #bin data into 4 bins reflecting likelihood of occupancy #create color palette plotcolors <- colorRampPalette(c("indianred3", "orange", "yellow", "green"))(4) #now, to put this on a google map map.CC.past <- get_map(location = c(lon = mean(df.CC.past$long), lat = mean(df.CC.past$lat)), zoom = 5, maptype = "terrain-background", source = "stamen" ) plotmap.CC.past<-ggmap(map.CC.past, extent="panel") plotmap.CC.past<-plotmap.CC.past+ geom_point(data=df.CC.past, aes(x=long, y=lat, color=value))+ annotate("text", x=-116.6, y=32.4, label="© Stamen Design, under CC BY 3.0. Data © OpenStreetMap, under ODbL.", size=3.1)+ scale_color_manual(values= plotcolors, name="Probability \nof Occupancy", labels=c("10 - 27%","27% - 44%","44 - 61%","61 - 78%")) plotmap.CC.past ggsave("CalCurrent_past_plot_ggmaps.pdf") myGLMs.past <- BIOMOD_LoadModels(myBiomodModelOut, models=c('GLM')) myRespPlot2D.past <- response.plot2(models = myGLMs.past, Data = get_formal_data(myBiomodModelOut, 'expl.var'), show.variables= get_formal_data(myBiomodModelOut,'expl.var.names'), do.bivariate = F, display_title= F, fixed.var.metric = 'mean', legend = F, data_species = get_formal_data(myBiomodModelOut, 'resp.var'))