#load libraries
library(raster)
library(dismo)
library(ncdf4)
library(dplyr)
library(matrixStats)
#https://stackoverflow.com/questions/17376939/problems-when-trying-to-load-a-package-in-r-due-to-rjava
Sys.setenv(JAVA_HOME = "C:/Program Files/Java/jre-1.8/")

setwd("D:/PhD/articles/article SDM/antarctic/BioOracle3/")


species <- c(
  "Alluroteuthis antarcticus",
  "Bathyteuthis abyssicola",
  "Galiteuthis glacialis",
  "Gonatus antarcticus",
  "Histioteuthis atlantica",
  "Histioteuthis eltaninae",
  "Kondakovia longimana",
  "Martialia hyadesi",
  "Mesonychoteuthis hamiltoni",
  "Moroteuthis ingens",
  "Moroteuthis robsoni",
  "Psychroteuthis glacialis",
  "Slosarczykovia circumantarctica",
  "Teuthowenia pellucida",
  "Todarodes filippovae"
)

for (species_name in species){
#species test
#species_name <- 'Alluroteuthis antarcticus'
load(file = paste(species_name,'Mdata.Rdata',sep=''))

  W = c(0,0,0,0,0)
  for (i in 1:5){
    if (AUCs[i]-.7>0){
      W[i]<- (100-100*abs(AUCs[i]-.9))^2
    }else {
      W[i]<-0
    }
  }
  W <- W/sum(W)
#calculate threshold
t1 <- threshold(em1,stat='spec_sens')
t2 <- threshold(em2,stat='spec_sens')
t3 <- threshold(em3,stat='spec_sens')
t4 <- threshold(em4,stat='spec_sens')
t5 <- threshold(em5,stat='spec_sens')
thresh <- t1*W[1]+t2*W[2]+t3*W[3]+t4*W[4]+t5*W[5]
#Analysis
#latitudinal Habitat SUitability
#line graphs, with +/-sd
rasterIn <- raster(paste(species_name,'_Present.tif',sep=''))
YY <-ymax(extent(rasterIn))

rasterIn <- as.matrix(rasterIn)
# r_m <- rowQuantiles(rasterIn,probs=.5,na.rm=TRUE)
r_m <- rowMeans(rasterIn,na.rm=TRUE)
r_s <- r_m+rowSds(rasterIn,na.rm=TRUE)
r_M <- apply(rasterIn, 1, max, na.rm=TRUE)
# row_25 <- rowQuantiles(rasterIn,probs=.25,na.rm=TRUE)#https://stackoverflow.com/questions/64132103/what-is-the-fastest-way-to-find-quantiles-of-each-row-in-a-matrix-with-r
# row_75 <- rowQuantiles(rasterIn,probs=.75,na.rm=TRUE)

svg(file=paste(species_name,"_lines_scenarios.svg",sep=''), width=10, height=10, onefile=TRUE,pointsize = 17)
par(mar = c(1, 1, 1, 2))
#present
plot(r_m, -.05*(1:length(r_m))+YY,
	type="l",col = "black", ylab = "Latitude" ,
	xlab = "Habitat suitability change", lwd=5.0, xlim=c(-0.1,1), axes=FALSE)
# lines(row_25, -.05*(1:length(r_m))+YY,col = "darkgrey", lwd=1.0, lty=2)
# lines(row_75, -.05*(1:length(r_m))+YY,col = "darkgrey", lwd=1.0, lty=2)
lines(r_s, -.05*(1:length(r_m))+YY,col = "black", lwd=3.0, lty=3)
lines(r_M, -.05*(1:length(r_m))+YY,col = "black", lwd=2.5, lty=3)

#2050
# #26
# rasterIn <- raster(paste("2050","26",paste(species_name,"2050","26",'.tif',sep=''),sep='/'))
# rasterIn <- as.matrix(rasterIn)
# r_m <- rowQuantiles(rasterIn,probs=.5,na.rm=TRUE)
# row_25 <- rowQuantiles(rasterIn,probs=.25,na.rm=TRUE)
# row_75 <- rowQuantiles(rasterIn,probs=.75,na.rm=TRUE)
# lines(r_m, -.05*(1:length(r_m))+YY,col = "blue", lwd=2.0, lty=1)
# lines(row_25, -.05*(1:length(r_m))+YY,col = "blue", lwd=1.0, lty=2)
# lines(row_75, -.05*(1:length(r_m))+YY,col = "blue", lwd=1.0, lty=2)
# 
# #85
rasterIn <- raster(paste("2050","85",paste(species_name,"2050","85",'.tif',sep=''),sep='/'))
rasterIn <- as.matrix(rasterIn)
r_m <- rowMeans(rasterIn,na.rm=TRUE)
r_s <- r_m+rowSds(rasterIn,na.rm=TRUE)
r_M <- apply(rasterIn, 1, max, na.rm=TRUE)
lines(r_m, -.05*(1:length(r_m))+YY,col = "purple", lwd=4.0)#, lty=5
# lines(r_m, -.05*(1:length(r_m))+YY,col = "red", lwd=2.0, lty=1)
# lines(row_25, -.05*(1:length(r_m))+YY,col = "red", lwd=1.0, lty=2)
# lines(row_75, -.05*(1:length(r_m))+YY,col = "red", lwd=1.0, lty=2)
lines(r_s, -.05*(1:length(r_m))+YY,col = "purple", lwd=2.5, lty=3)
lines(r_M, -.05*(1:length(r_m))+YY,col = "purple", lwd=1.5, lty=3)

#2100
# #26
# rasterIn <- raster(paste("2100","26",paste(species_name,"2100","26",'.tif',sep=''),sep='/'))
# rasterIn <- as.matrix(rasterIn)
# r_m <- rowMeans(rasterIn,na.rm=TRUE)
# r_s <- r_m+rowSds(rasterIn,na.rm=TRUE)
# r_M <- apply(rasterIn, 1, max, na.rm=TRUE)
# # r_m <- rowQuantiles(rasterIn,probs=.5,na.rm=TRUE)
# # row_25 <- rowQuantiles(rasterIn,probs=.25,na.rm=TRUE)
# # row_75 <- rowQuantiles(rasterIn,probs=.75,na.rm=TRUE)
# lines(r_m, -.05*(1:length(r_m))+YY,col = "blue", lwd=1.5)#, lty=5
# # lines(row_25, -.05*(1:length(r_m))+YY,col = "blue", lwd=1.0, lty=3)
# # lines(row_75, -.05*(1:length(r_m))+YY,col = "blue", lwd=1.0, lty=3)
# lines(r_s, -.05*(1:length(r_m))+YY,col = "blue", lwd=1.0, lty=3)
# lines(r_M, -.05*(1:length(r_m))+YY,col = "blue", lwd=1.0, lty=3)

#85
rasterIn <- raster(paste("2100","85",paste(species_name,"2100","85",'.tif',sep=''),sep='/'))
rasterIn <- as.matrix(rasterIn)
r_m <- rowMeans(rasterIn,na.rm=TRUE)
r_s <- r_m+rowSds(rasterIn,na.rm=TRUE)
r_M <- apply(rasterIn, 1, max, na.rm=TRUE)
# r_m <- rowQuantiles(rasterIn,probs=.5,na.rm=TRUE)
# row_25 <- rowQuantiles(rasterIn,probs=.25,na.rm=TRUE)
# row_75 <- rowQuantiles(rasterIn,probs=.75,na.rm=TRUE)
lines(r_m, -.05*(1:length(r_m))+YY,col = "green", lwd=4.0)#, lty=5
# lines(row_25, -.05*(1:length(r_m))+YY,col = "red", lwd=1.0, lty=3)
# lines(row_75, -.05*(1:length(r_m))+YY,col = "red", lwd=1.0, lty=3)
lines(r_s, -.05*(1:length(r_m))+YY,col = "green", lwd=2.5, lty=3)
lines(r_M, -.05*(1:length(r_m))+YY,col = "green", lwd=1.5, lty=3)

segments(thresh, -30, thresh, -79, lwd=2.0, col = "darkgrey")
axis(1)   # Draw the x-axis below the plot area
axis(4)   # Draw the y-axis to the right of the plot area
box()
dev.off()


#Change of total Habitat suitability
#Habitat suitability changed to threshold map

sink(paste(species_name, "_area.txt", sep=''))
cat(paste(species_name,'\n'))
#apply to suitability maps
rasterIn <- raster(paste(species_name,'_Present.tif',sep=''))
A <- area(rasterIn)
map <- rasterIn>thresh
writeRaster(map,file=paste(species_name,'_Present_binary.tif',sep=''), options="INTERLEAVE=BAND", overwrite=TRUE)

#total area of habitat
cat("present\n")
cat(paste(sum(as.matrix(A)*as.matrix(map),na.rm=TRUE),'\n'))

for (RCP in c('26','85')){
	for (time in c('2050','2100')){
  # time <- '2050'
	  cat(paste(time,RCP,'\n'))
		#read raster
		rasterIn <- raster(paste(time,RCP,paste(species_name,time,RCP,'.tif',sep=''),sep='/'))
		#apply threshold
		map <- rasterIn>thresh
		#save raster with threshold for hotspot map
		writeRaster(map,file=paste(time,RCP,paste(species_name,time,RCP,'_binary.tif',sep=''),sep='/') , options="INTERLEAVE=BAND", overwrite=TRUE)
    writeRaster(map-raster(paste(species_name,'_Present_binary.tif',sep='')),file=paste(time,RCP,paste(species_name,time,RCP,'_binary_diff.tif',sep=''),sep='/') , options="INTERLEAVE=BAND", overwrite=TRUE)
    #total area of habitat
		cat(paste(sum(as.matrix(A)*as.matrix(map),na.rm=TRUE),'\n'))
	 }
}
sink()

}





#Sum all species for biodiversity Hotspot
for (RCP in c('26','85')){
  for (time in c('2050','2100')){
    rasterIn <- raster(paste(time,RCP,paste("Alluroteuthis antarcticus",time,RCP,'_binary.tif',sep=''),sep='/'))
    values(rasterIn) <- 0
    for (species_name in species){
      rasterIn <- rasterIn+raster(paste(time,RCP,paste(species_name,time,RCP,'_binary.tif',sep=''),sep='/'))
    }
    writeRaster(rasterIn,file=paste(time,RCP,paste('biodiversity',time,RCP,'.tif',sep=''),sep='/') , options="INTERLEAVE=BAND", overwrite=TRUE)
  }
}
#Present
rasterIn <- raster(paste("Alluroteuthis antarcticus",'_Present_binary.tif',sep=''))
values(rasterIn) <- 0
for (species_name in species){
  rasterIn <- rasterIn+raster(paste(species_name,'_Present_binary.tif',sep=''))
}
writeRaster(rasterIn,file=paste('biodiversity_Present_binary.tif',sep=''), options="INTERLEAVE=BAND", overwrite=TRUE)




#Sum all diffs but filter for negative and positive diffs
for (RCP in c('26','85')){
  for (time in c('2050','2100')){
    mapi <- raster(paste(time,RCP,paste("Alluroteuthis antarcticus",time,RCP,'_binary.tif',sep=''),sep='/'))
    values(mapi) <- 0
    mapii <- mapi
    #duplicate for 2 kinds of differences (positive and negative)
    for (species_name in species){
      rasterIn <-raster(paste(time,RCP,paste(species_name,time,RCP,'_binary_diff.tif',sep=''),sep='/'))
      #negative
      #writeRaster
      mapi <- mapi+(rasterIn<0)
      #positive
      #writeRaster
      mapii <- mapii+(rasterIn>0)
    }
    writeRaster(mapi,file=paste(time,RCP,paste('negative_species.tif',sep=''),sep='/'), options="INTERLEAVE=BAND", overwrite=TRUE)
    writeRaster(mapii,file=paste(time,RCP,paste('positive_species.tif',sep=''),sep='/'), options="INTERLEAVE=BAND", overwrite=TRUE)
  }
}