## written by Olivier Broennimann. Departement of Ecology and Evolution (DEE). ## 16.01.2019. University of Lausanne. Department of ecology and evolution, Switzerland t1=Sys.time() #install.packages("devtools") #library(devtools) #install_github(repo="ecospat/ecospat/ecospat",ref="3.1") library(ecospat) library(ade4) library(raster) library(rworldmap) ########################### #### data preparation ##### ########################### # get climate data for the two area clim<-getData('worldclim', var='bio', res=10) # choice of variables var<-c("bio1","bio5","bio7","bio12","bio16","bio18") clim<-clim[[which(names(clim)%in%var)]] # create mask for native area background data(countriesCoarseLessIslands) #countries from rworldmap ctry = c("AUT","BEL","BGR","HRV","CYP","CZE","DNK","EST","FIN","FRA","DEU","GRC","HUN","IRL", "ITA","LVA","LTU","LUX","MLT","NLD","POL","PRT","ROU","SVK","SVN","ESP","SWE","GBR", "NOR","CHE","BLR","UKR","MDA","MNE","SRB","BIH","ALB","MKD") bkg.eu<-aggregate(countriesCoarseLessIslands[countriesCoarseLessIslands$ADM0_A3%in%ctry,]) # create mask for invaded area background ctry = c("CAN", "USA") bkg.nam<-aggregate(countriesCoarseLessIslands[countriesCoarseLessIslands$ADM0_A3%in%ctry,]) # extract climate data from the rasters clim.bkg.eu<-mask(crop(clim,bbox(bkg.eu)),bkg.eu) clim.bkg.nam<-mask(crop(clim,bbox(bkg.nam)),bkg.nam) # load occurrence sites for the species (column names should be x,y) file<-"https://www.unil.ch/ecospat/files/live/sites/ecospat/files/shared/Files_site/CStoebe_EU_NA.txt" occ.sp.aggr.dupl<-read.delim(file,h=T,sep="\t") # remove duplicated occurrences occ.sp.aggr<-occ.sp.aggr.dupl[!duplicated(occ.sp.aggr.dupl),] # remove occurrences closer than a minimum distance to each other (remove aggregation). Setting min.dist=0 will remove no occurrence. If your dataset is in degree, 0.008333 correspond to 1km at the equator. occ.sp<-ecospat.occ.desaggregation(occ.sp.aggr,min.dist = 1/60*10) occ.eu<-base::subset(occ.sp,region=="EU",select = c(x,y)) occ.nam<-base::subset(occ.sp,region=="NAM",select = c(x,y)) # create sp occurrence dataset by extracting climate variables from the rasters env.bkg.eu<-na.exclude(data.frame(extract(clim,bkg.eu))) env.bkg.nam<-na.exclude(data.frame(extract(clim,bkg.nam))) env.occ.eu<-na.exclude(extract(clim,occ.eu)) env.occ.nam<-na.exclude(extract(clim,occ.nam)) #check spatial autocorrelation ecospat.mantel.correlogram(dfvar=cbind(occ.eu,env.occ.eu),colxy=1:2,n=200,nclass=10,nperm=20) ecospat.mantel.correlogram(dfvar=cbind(occ.nam,env.occ.nam),colxy=1:2,n=200,nclass=10,nperm=20) ################################ #### niche quantifications ##### ################################ #calibration of PCA-env pca.env <-dudi.pca(rbind(env.bkg.eu,env.bkg.nam), center = T, scale = T, scannf = F, nf = 2) # predict the scores on the PCA axes scores.bkg<- pca.env$li scores.bkg.eu<- suprow(pca.env,env.bkg.eu)$lisup scores.bkg.nam<- suprow(pca.env,env.bkg.nam)$lisup scores.occ.eu<- suprow(pca.env,env.occ.eu)$lisup scores.occ.nam<- suprow(pca.env,env.occ.nam)$lisup # calculation of occurence density z1<- ecospat.grid.clim.dyn(scores.bkg,scores.bkg.eu,scores.occ.eu,R=100) z2<- ecospat.grid.clim.dyn(scores.bkg,scores.bkg.nam,scores.occ.nam,R=100) equ<-ecospat.niche.equivalency.test(z1,z2,rep=10) # test of niche equivalency and similarity sim1<-ecospat.niche.similarity.test(z1,z2,rep=100,alternative = "greater",rand.type = 1) #niches randomly shifted in both area sim2<-ecospat.niche.similarity.test(z1,z2,rep=100,alternative = "greater",rand.type = 2) #niche randomly shifted only in invaded area # overlap corrected by availabilty of background conditions ecospat.niche.overlap(z1,z2,cor=T) # uncorrected overlap ecospat.niche.overlap(z1,z2,cor=F) ################################ #### niche visualizations ##### ################################ # occurrence density plots ecospat.plot.niche(z1,title="PCA-env - EU niche",name.axis1="PC1",name.axis2="PC2") ecospat.plot.niche(z2,title="PCA-env - NA niche",name.axis1="PC1",name.axis2="PC2") # contribution of original variables ecospat.plot.contrib(pca.env$co,pca.env$eig) # niche tests plots ecospat.plot.overlap.test(equ,"D","Equivalency") ecospat.plot.overlap.test(sim1,"D","Similarity 1<->2") #niches randomly shifted in both areas ecospat.plot.overlap.test(sim2,"D","Similarity 1->2") #niche randomly shifted only in invaded area ecospat.plot.niche.dyn(z1,z2,quant=0.5,name.axis1="PC1",name.axis2="PC2",interest=1) ecospat.plot.niche.dyn(z1,z2,quant=0.5,name.axis1="PC1",name.axis2="PC2",interest=2) # occurrence density in geography plot(ecospat.niche.zProjGeo(z1,clim.bkg.eu)) points(occ.eu,cex=0.2,pch=19) plot(ecospat.niche.zProjGeo(z2,clim.bkg.nam)) points(occ.nam,cex=0.2,pch=19) t2=Sys.time() t2-t1 #Time difference of 1.554178 mins