library(SDMTools)#load necessary libraries
future.dir= 'E:/SEABat/new/future/'
current.dir = 'E:/SEABat/new/current/'; setwd(current.dir) #define & set the working directory - current
out.dir = 'E:/SEABat/new/overlap/';dir.create(out.dir);

###01.create a vector of species names
files=list.files(pattern='_binary') #get a list of files in your working directory that include 'asc'
files=grep('.asc',files,value=T)

species=NULL #create the object name of your vector
spp=NULL #create a temporary object name for elements of your vector

for(file in files) { cat(file,'\n') #loop through the files, and for each file, complete the following:
spp=strsplit(file,'\\.') #split your file name by '.'  the two resulting elements will be 'Rhinolophus_sp' and 'asc'
spp=spp[[1]][1] #select the first of your two elements and call it 'spp'.  spp will be 'Rhinolophus_sp'

if(length(species)==0){ #if species is still null, ie no species in your vector yet
	species=spp #then add the first species
	} else { #once you have one species already in the vector
	species=c(species,spp)} #bind each new species to it
}

cols = c('gray86','red3','forestgreen','blue')
overlap=NULL
overlap.matrix=NULL
out=NULL
spp_name=NULL
for (spp in species) {
#setup some plot parameters
legend.pnts = cbind(c(93,91,91,93),c(2,2,7.5,7.5))

#read in the current binary ascii
current = read.asc.gz(paste(spp,'.asc.gz',sep=''))

cs.current=ClassStat(current, latlon=TRUE)
cs.current$total.area = cs.current$total.area / 1000000000 #area, 1000s km^2
	  
#read in the future binary ascii
future = read.asc.gz(paste(future.dir, spp,'.asc.gz',sep=''))
future[which(is.finite(future) & future>threshold)] = 2 #set all parts of the species distribuiton to 2
future[which(future<1)]=0

cs.future=ClassStat(future, latlon=TRUE)
cs.future$total.area = cs.future$total.area / 1000000000 #area, 1000s km^2

overlap=current+future

cs.overlap=ClassStat(overlap, latlon=TRUE)
cs.overlap$total.area = cs.overlap$total.area / 1000000000 #area, 1000s km^2

out=matrix(NA,nr=1,ncol=3)
out=as.data.frame(out)
colnames(out)=c('current','future','overlap')
out$current=cs.current[2,3]
out$future=cs.future[2,3]
out$overlap=cs.overlap[4,3]

write.csv(out, paste(out.dir, spp, '_area.csv'))
write.asc.gz(overlap,paste(out.dir, spp, '_overlap.asc', sep=''))

spp_name=strsplit(spp,'_')[[1]][1:2]
spp_name=paste(spp_name[1], '_', spp_name[2])

if (is.null(overlap.matrix)){
	overlap.matrix=cbind(spp_name,out)
	} else {
	overlap.matrix=rbind(overlap.matrix,cbind(spp_name,out))
	}
	

png(paste(out.dir, spp,'_overlap.png', sep=''), width=nrow(bat.asc),height=ncol(bat.asc), units='px', res=300, pointsize=5, bg='white') #start the plot
	par(mar=c(0,0,0,0)+0.1) #remove any plot margins
	image(overlap, ann=FALSE,axes=FALSE,col=cols, zlim=c(0,3)) #plot the richness
	legend.gradient(legend.pnts,cols=cols,limits=c(0,3), title='overlap', cex=2)
dev.off()
}

write.csv(overlap.matrix, paste(out.dir, "all_species.csv",sep=''))