#Step 3 - Calculate means of MAXENT replicates #This script calculates means of replicate MAXENT models and saves them in a separate folder. #This takes some time. library(raster) library(dismo) library(rgdal) library(parallel) setwd("") #Set the directory where MAXENT projections are stored (during Step 2) st.dir <- "/maxentRpred" st.bigdir <- "/maxentRpredmeans" #If you're having trouble loading the rasters... rasterOptions(tmpdir = "G:/Rtemp/") #List the filenames of the MAXENT projections ch.files <- list.files(st.dir) #Keep only .tif filenames ch.files <- ch.files[sapply(ch.files, function(i) substr(i, nchar(i) - 3, nchar(i)) == ".tif")] #Get positions of hyphens, in order to extract metadata encoded in filename mx.pos <- t(sapply(ch.files, function(i) c(as.vector(gregexpr(pattern = "-", i)[[1]]), nchar(i)))) #Translate file metadata (e.g. species id, region) into a matrix mx.files <- t(sapply(1:length(ch.files), function(i) c( ch.files[i], substr(ch.files[i], mx.pos[i, 1] + 1, mx.pos[i, 2] - 1), substr(ch.files[i], mx.pos[i, 2] + 1, mx.pos[i, 3] - 1), substr(ch.files[i], mx.pos[i, 3] + 1, mx.pos[i, 4] - 1), substr(ch.files[i], mx.pos[i, 5] + 1, mx.pos[i, 6] - 1), substr(ch.files[i], mx.pos[i, 6] + 1, mx.pos[i, 7] - 4) ))) #Remove duplicates mx.noreps <- mx.files[duplicated(mx.files[,2:5]) == F,] #The following steps help to only analyze models that have not already had means calculated ch.need <- sapply(1:nrow(mx.noreps), function(i) { v.i <- mx.noreps[i,] paste("pred-", v.i[2], "-", v.i[3], "-", v.i[4], "-", v.i[5], ".tif", sep = "") }) ch.done <- list.files(st.bigdir) v.done <- unname(sapply(ch.done, function(i) which(ch.need == i))) #Only calculate means of projections that have not already had means calculated v.do <- 1:nrow(mx.noreps) v.do <- v.do[-v.done] #This will calculate means of the 10 replicate projections of each model, and save in the maxentRpredmeans folder s.time <- proc.time() cl <- makeCluster(6) clusterEvalQ(cl, {library(dismo); library(rgdal); library(raster)}) clusterExport(cl, c("mx.files", "mx.noreps", "st.dir", "v.do")) l.done <- parLapply(cl, v.do, function(i) { rasterOptions(tmpdir = "G:/Rtemp/") v.i <- mx.noreps[i,] v.use <- mx.files[ mx.files[,2] == v.i[2] & mx.files[,3] == v.i[3] & mx.files[,4] == v.i[4] & mx.files[,5] == v.i[5], 1 ] r.use <- mean(stack(sapply(v.use, function(j) raster(paste(st.dir, "/", j, sep = ""))))) writeRaster( r.use, paste(st.bigdir, v.i[2], "-", v.i[3], "-", v.i[4], "-", v.i[5], ".tif", sep = ""), overwrite = T ) rm(r.use) gc() i }) stopCluster(cl) proc.time() - s.time