library(raster) library (maptools) # load in lat/lon data setwd ("../Dropbox/squamate_speciation/R") DATA <- read.csv ("dataset.csv") SPECIES <- DATA[which(DATA$Pair==1),] SPECIES <- transform (SPECIES, ID=as.numeric(factor(species))) LONGLAT <- subset (SPECIES)[,c("longitude","latitude")] #get worldclim data from internet if you want 10, 5, or 2.5 min. If you need 30 sec, you need to provide tile coordinates lon=5, lat=45 #climatelayers = getData('worldclim', var='bio', res=2.5, path=tempdir()) VARIABLES <- list.files (path="C:/Users/Tereza/Dropbox/wordclim/bio_30s", pattern='bil', full.names=TRUE) WC30 <- stack (VARIABLES) #background points SPECIES1 <- SPECIES[which(SPECIES$ID==1),] SPECIES1 <- SPECIES1 [,3:4] SPECIES2 <- SPECIES[which(SPECIES$ID==2),] SPECIES2 <- SPECIES2 [,3:4] PAIR <- SPECIES [,3:4] library (dismo) hull <- convHull(PAIR, lonlat=TRUE) hull1 <- convHull(SPECIES1, lonlat=TRUE) hull2 <- convHull(SPECIES2, lonlat=TRUE) hullbg <- hull@polygons - (hull1@polygons + hull2@polygons) crs (hullbg) <- crs (wrld_simpl) CONTINENT <- readShapePoly ("continents/continent.shp") crs (CONTINENT) <- crs (wrld_simpl) hullbg <- intersect (hullbg, CONTINENT) library (rgeos) hullbg <- gBuffer (hullbg, width=-0.1) crs (hullbg) <- crs (wrld_simpl) data (wrld_simpl) plot (hull@polygons, border='yellow', lwd=2) points (SPECIES1, cex=1, col='red') points (SPECIES2, cex=1, col='blue') plot (hull1@polygons, border='red', lwd=2, add=TRUE) plot (hull2@polygons, border='blue', lwd=2, add=TRUE) plot (wrld_simpl, add=TRUE) plot (hullbg, border='green', lwd=2, add=TRUE) NRECORDS <- nrow (SPECIES) box () # extract transformed climate values CLIMATE <- extract(WC30, LONGLAT) CLIMATE_RECORDS <- cbind (SPECIES, CLIMATE) CLIMATE_RECORDS <- (CLIMATE_RECORDS [, 5:24]) CLIMATE_RECORDS <- CLIMATE_RECORDS[complete.cases(CLIMATE_RECORDS), ] ID <- CLIMATE_RECORDS [,1] CLIMATE <- CLIMATE_RECORDS [,-1] replicate (100, { BG <- spsample (hullbg, n=NRECORDS, "random") BG_CLIMATE <- extract(WC30, BG) ALL_CLIMATE <- rbind (CLIMATE, BG_CLIMATE) BG_ID <- transform (BG_CLIMATE, ID=3) BG_ID <- BG_ID [, 20] ID <- c(ID, BG_ID) NRECORDS1 <- nrow (SPECIES1) NRECORDS1 NRECORDS2 <- nrow (SPECIES2) NRECORDS2 PCA <- prcomp (ALL_CLIMATE, center=TRUE, scale.=TRUE) library (vegan) EIGEN<-eigenvals (PCA) if (EIGEN["PC2"]<1) { NUM <- 1 } else if (EIGEN["PC3"]<1) { NUM <- 2 } else if (EIGEN["PC4"]<1) { NUM <- 3 } else if (EIGEN["PC5"]<1) { NUM <- 4 } else { NUM <- 5 } if (NRECORDS1(NUM2/10)) { NUM_FINAL <- (round (NUM2/10, digits=0) +1) } else { NUM_FINAL <- NUM+1 } #plot(PCA, type = "l") SCORES <- PCA$x SCORES <- cbind (ID, SCORES) SCORES1 <- subset(SCORES, ID=="1") SCORES1 <- SCORES1 [,2:3] SCORES2 <- subset(SCORES, ID=="2") SCORES2 <- SCORES2 [,2:3] SCORES3 <- subset(SCORES, ID=="3") SCORES3 <- SCORES3 [,2:3] # compute hypervolumes with auto-bandwidth for both species library (hypervolume) hv1 <- hypervolume(SCORES1,quantile=0.1,reps=1000,bandwidth=estimate_bandwidth(SCORES1),name='1') hv2 <- hypervolume(SCORES2,quantile=0.1,reps=1000,bandwidth=estimate_bandwidth(SCORES2),name='2') hv3 <- hypervolume(SCORES3,quantile=0.1,reps=1000,bandwidth=estimate_bandwidth(SCORES3),name='bg') # determine intersection and unique components of the overlap hv_set12 <- hypervolume_set(hv1, hv2, check_memory=FALSE) hv_set13 <- hypervolume_set(hv1, hv3, check_memory=FALSE) hv_set23 <- hypervolume_set(hv2, hv3, check_memory=FALSE) # put all the output volumes in one convenient place volumes12 <- get_volume(hv_set12) volumes13 <- get_volume(hv_set13) volumes23 <- get_volume(hv_set23) PAIR <- SPECIES$Pair [1] volumes <- c (PAIR, NRECORDS1, NRECORDS2, NUM_FINAL, volumes12 ["HV1"], volumes12 ["HV2"], volumes12 ["Intersection"],volumes13 ["HV1"], volumes13 ["HV2"], volumes13 ["Intersection"], volumes23 ["HV1"], volumes23 ["HV2"], volumes23 ["Intersection"] ) round (volumes, digits=0) #rownames (volumes) <- c ("pair", "records_species1", "records_species2", "PCs","HV1", "HV2", "Intersection12", "HV1b", "HV3", "Intersection13", "HV2b", "HV3b", "Intersection23") volumes <- t (volumes) volumes #write.table(volumes, file = "volumes.csv", append=FALSE, sep = ",", row.names = FALSE,col.names = TRUE) write.table(volumes, file = "volumes.csv", append=TRUE, sep = ",", row.names = FALSE,col.names = FALSE) VOLUMES <- read.csv ("volumes.csv") VOLUMES_MEAN <- aggregate (.~pair, VOLUMES, function (x) c (mean = mean (x))) write.table(VOLUMES_MEAN, file = "volumes_mean.csv", append=FALSE, sep = ",", row.names = FALSE,col.names = TRUE) }) VOLUMES_MEAN <- read.csv ("volumes_mean.csv") VOLUMES_MEAN