###Latitudinal trends and changes in habitat suitability

#1. Import habitat suitability ensembles 

#Present
library(raster)

present <- raster("present_ensembled.bil")
projection(present)

#Future
f2050r2.6 <- raster("2050_RCP2.6_ensembled.bil")
f2050r4.5 <- raster("2050_RCP4.5_ensembled.bil")
f2050r6.0 <- raster("2050_RCP6.0_ensembled.bil")
f2050r8.5 <- raster("2050_RCP8.5_ensembled.bil")
f2100r2.6 <- raster("2100_RCP2.6_ensembled.bil")
f2100r4.5 <- raster("2100_RCP4.5_ensembled.bil")
f2100r6.0 <- raster("2100_RCP6.0_ensembled.bil")
f2100r8.5 <- raster("2100_RCP8.5_ensembled.bil")


#####################latitudinal mean habitat suitability####

#1st. manipulate the raster to obtain the yvals and xvals for each raster of interest

#Present
stepsize_present = (present@extent@ymax-present@extent@ymin)/present@nrows
yvals_present = seq(present@extent@ymax-stepsize_present / 2, present@extent@ymin, -stepsize_present)
xvals_present = rowMeans(as.matrix(present),na.rm=TRUE)

#205026
stepsize_205026 = (f2050r2.6@extent@ymax-f2050r2.6@extent@ymin)/f2050r2.6@nrows
yvals_205026 = seq(f2050r2.6@extent@ymax-stepsize_205026 / 2, f2050r2.6@extent@ymin, -stepsize_205026)
xvals_205026 = rowMeans(as.matrix(f2050r2.6),na.rm=TRUE)

#210026
stepsize_210026 = (f2100r2.6@extent@ymax-f2100r2.6@extent@ymin)/f2100r2.6@nrows
yvals_210026 = seq(f2100r2.6@extent@ymax-stepsize_210026 / 2, f2100r2.6@extent@ymin, -stepsize_210026)
xvals_210026 = rowMeans(as.matrix(f2100r2.6),na.rm=TRUE)

#205045
stepsize_205045 = (f2050r4.5@extent@ymax-f2050r4.5@extent@ymin)/f2050r4.5@nrows
yvals_205045 = seq(f2050r4.5@extent@ymax-stepsize_205045 / 2, f2050r4.5@extent@ymin, -stepsize_205045)
xvals_205045 = rowMeans(as.matrix(f2050r4.5),na.rm=TRUE)

#210045
stepsize_210045 = (f2100r4.5@extent@ymax-f2100r4.5@extent@ymin)/f2100r4.5@nrows
yvals_210045 = seq(f2100r4.5@extent@ymax-stepsize_210045 / 2, f2100r4.5@extent@ymin, -stepsize_210045)
xvals_210045 = rowMeans(as.matrix(f2100r4.5),na.rm=TRUE)

#205060
stepsize_205060 = (f2050r6.0@extent@ymax-f2050r6.0@extent@ymin)/f2050r6.0@nrows
yvals_205060 = seq(f2050r6.0@extent@ymax-stepsize_205060 / 2, f2050r6.0@extent@ymin, -stepsize_205060)
xvals_205060 = rowMeans(as.matrix(f2050r6.0),na.rm=TRUE)

#210060
stepsize_210060 = (f2100r6.0@extent@ymax-f2100r6.0@extent@ymin)/f2100r6.0@nrows
yvals_210060 = seq(f2100r6.0@extent@ymax-stepsize_210060 / 2, f2100r6.0@extent@ymin, -stepsize_210060)
xvals_210060 = rowMeans(as.matrix(f2100r6.0),na.rm=TRUE)

#205085
stepsize_205085 = (f2050r8.5@extent@ymax-f2050r8.5@extent@ymin)/f2050r8.5@nrows
yvals_205085 = seq(f2050r8.5@extent@ymax-stepsize_205085 / 2, f2050r8.5@extent@ymin, -stepsize_205085)
xvals_205085 = rowMeans(as.matrix(f2050r8.5),na.rm=TRUE)

#210085
stepsize_210085 = (f2100r8.5@extent@ymax-f2100r8.5@extent@ymin)/f2100r8.5@nrows
yvals_210085 = seq(f2100r8.5@extent@ymax-stepsize_210085 / 2, f2100r8.5@extent@ymin, -stepsize_210085)
xvals_210085 = rowMeans(as.matrix(f2100r8.5),na.rm=TRUE)

#2nd. Subtract the 2100/2050 matrices to the present-day matrix in each scenario:

#what we want is to subtract the x vals (the hab suit date). 
#Let's try it:

xvals_dif205026 <- xvals_205026 - xvals_present
xvals_dif205045 <- xvals_205045 - xvals_present
xvals_dif205060 <- xvals_205060 - xvals_present
xvals_dif205085 <- xvals_205085 - xvals_present

xvals_dif210026 <- xvals_210026 - xvals_present
xvals_dif210045 <- xvals_210045 - xvals_present
xvals_dif210060 <- xvals_210060 - xvals_present
xvals_dif210085 <- xvals_210085 - xvals_present

ppi=360

#205026
x.matrix205026 = matrix(data = xvals_dif205026, ncol = 1, byrow = TRUE)
y.matrix = matrix(data = yvals_present, ncol = 1, byrow = TRUE)
df_205026 <- cbind(x.matrix205026,y.matrix)
class(try)
df_205026 <- as.data.frame(df_205026)
library(dplyr)
df_205026 <- df_205026 %>% 
  mutate(C = if_else(V1 > 0, "Gain", "Loss"))
df_205026$C <- as.factor(df_205026$C)

library(ggplot2)
plot205026 <- ggplot(data=df_205026, aes(x=V2, y=V1, fill = C)) +
  geom_bar(stat="identity")+
  geom_hline(yintercept = 0, linetype="dashed", 
             color = "black", size=1.2)+
  ylim(-0.022,0.022)+
  scale_fill_manual(values=c("red4",
                             "royalblue4"))+
  labs(x = "Latitude (meters)", y = "Difference in mean habitat suitability (%)")+
  theme(axis.title.y = element_text(angle=90, vjust=2.125, face="bold", size=18),
        axis.text.y = element_text(size=14)) +
  theme(panel.grid = element_line(color = "darkgrey",
                                  size = 0.3,
                                  linetype = 2))+
  scale_x_continuous(n.breaks = 20)+
  theme(axis.text.x = element_text(size=14),axis.title.x = element_text(angle=0, vjust = 0, face="bold", size=18))+
  
  theme(panel.border = element_rect(size=0.7, colour = "black", fill=NA, linetype=1))+
  theme(legend.position="none")+
  ggtitle("P. fraudulenta, 2050_RCP-2.6")+
  coord_flip()

#210026
x.matrix210026 = matrix(data = xvals_dif210026, ncol = 1, byrow = TRUE)
y.matrix = matrix(data = yvals_present, ncol = 1, byrow = TRUE)
df_210026 <- cbind(x.matrix210026,y.matrix)
class(try)
df_210026 <- as.data.frame(df_210026)
library(dplyr)
df_210026 <- df_210026 %>% 
  mutate(C = if_else(V1 > 0, "Gain", "Loss"))
df_210026$C <- as.factor(df_210026$C)

library(ggplot2)
plot210026 <- ggplot(data=df_210026, aes(x=V2, y=V1, fill = C)) +
  geom_bar(stat="identity")+
  geom_hline(yintercept = 0, linetype="dashed", 
             color = "black", size=1.2)+
  ylim(-0.022,0.022)+
  scale_fill_manual(values=c("red4",
                             "royalblue4"))+
  labs(x = "Latitude (meters)", y = "Difference in mean habitat suitability (%)")+
  theme(axis.title.y = element_text(angle=90, vjust=2.125, face="bold", size=18),
        axis.text.y = element_text(size=14)) +
  theme(panel.grid = element_line(color = "darkgrey",
                                  size = 0.3,
                                  linetype = 2))+
  scale_x_continuous(n.breaks = 20)+
  theme(axis.text.x = element_text(size=14),axis.title.x = element_text(angle=0, vjust = 0, face="bold", size=18))+
  
  theme(panel.border = element_rect(size=0.7, colour = "black", fill=NA, linetype=1))+
  theme(legend.position="none")+
  ggtitle("P. fraudulenta, 2100_RCP-2.6")+
  coord_flip()


tiff(filename= 'fraudulenta_habdif205026.tiff',width = 12*ppi, height = 8*ppi, res = ppi)
plot205026
dev.off()
tiff(filename= 'fraudulenta_habdif210026.tiff',width = 12*ppi, height = 8*ppi, res = ppi)
plot210026
dev.off()

#205045
x.matrix205045 = matrix(data = xvals_dif205045, ncol = 1, byrow = TRUE)
y.matrix = matrix(data = yvals_present, ncol = 1, byrow = TRUE)
df_205045 <- cbind(x.matrix205045,y.matrix)
class(try)
df_205045 <- as.data.frame(df_205045)
library(dplyr)
df_205045 <- df_205045 %>% 
  mutate(C = if_else(V1 > 0, "Gain", "Loss"))
df_205045$C <- as.factor(df_205045$C)

library(ggplot2)
plot205045 <- ggplot(data=df_205045, aes(x=V2, y=V1, fill = C)) +
  geom_bar(stat="identity")+
  geom_hline(yintercept = 0, linetype="dashed", 
             color = "black", size=1.2)+
  ylim(-0.022,0.022)+
  scale_fill_manual(values=c("red4",
                             "royalblue4"))+
  labs(x = "Latitude (meters)", y = "Difference in mean habitat suitability (%)")+
  theme(axis.title.y = element_text(angle=90, vjust=2.125, face="bold", size=18),
        axis.text.y = element_text(size=14)) +
  theme(panel.grid = element_line(color = "darkgrey",
                                  size = 0.3,
                                  linetype = 2))+
  scale_x_continuous(n.breaks = 20)+
  theme(axis.text.x = element_text(size=14),axis.title.x = element_text(angle=0, vjust = 0, face="bold", size=18))+
  
  theme(panel.border = element_rect(size=0.7, colour = "black", fill=NA, linetype=1))+
  theme(legend.position="none")+
  ggtitle("P. fraudulenta, 2050_RCP-4.5")+
  coord_flip()

#210045
x.matrix210045 = matrix(data = xvals_dif210045, ncol = 1, byrow = TRUE)
y.matrix = matrix(data = yvals_present, ncol = 1, byrow = TRUE)
df_210045 <- cbind(x.matrix210045,y.matrix)
class(try)
df_210045 <- as.data.frame(df_210045)
library(dplyr)
df_210045 <- df_210045 %>% 
  mutate(C = if_else(V1 > 0, "Gain", "Loss"))
df_210045$C <- as.factor(df_210045$C)

library(ggplot2)
plot210045 <- ggplot(data=df_210045, aes(x=V2, y=V1, fill = C)) +
  geom_bar(stat="identity")+
  geom_hline(yintercept = 0, linetype="dashed", 
             color = "black", size=1.2)+
  ylim(-0.022,0.022)+
  scale_fill_manual(values=c("red4",
                             "royalblue4"))+
  labs(x = "Latitude (meters)", y = "Difference in mean habitat suitability (%)")+
  theme(axis.title.y = element_text(angle=90, vjust=2.125, face="bold", size=18),
        axis.text.y = element_text(size=14)) +
  theme(panel.grid = element_line(color = "darkgrey",
                                  size = 0.3,
                                  linetype = 2))+
  scale_x_continuous(n.breaks = 20)+
  theme(axis.text.x = element_text(size=14),axis.title.x = element_text(angle=0, vjust = 0, face="bold", size=18))+
  
  theme(panel.border = element_rect(size=0.7, colour = "black", fill=NA, linetype=1))+
  theme(legend.position="none")+
  ggtitle("P. fraudulenta, 2100_RCP-4.5")+
  coord_flip()

tiff(filename= 'fraudulenta_habdif205045.tiff',width = 12*ppi, height = 8*ppi, res = ppi)
plot205045
dev.off()
tiff(filename= 'fraudulenta_habdif210045.tiff',width = 12*ppi, height = 8*ppi, res = ppi)
plot210045
dev.off()

#205060
x.matrix205060 = matrix(data = xvals_dif205060, ncol = 1, byrow = TRUE)
y.matrix = matrix(data = yvals_present, ncol = 1, byrow = TRUE)
df_205060 <- cbind(x.matrix205060,y.matrix)
class(try)
df_205060 <- as.data.frame(df_205060)
library(dplyr)
df_205060 <- df_205060 %>% 
  mutate(C = if_else(V1 > 0, "Gain", "Loss"))
df_205060$C <- as.factor(df_205060$C)

library(ggplot2)
plot205060 <- ggplot(data=df_205060, aes(x=V2, y=V1, fill = C)) +
  geom_bar(stat="identity")+
  geom_hline(yintercept = 0, linetype="dashed", 
             color = "black", size=1.2)+
  ylim(-0.022,0.022)+
  scale_fill_manual(values=c("red4",
                             "royalblue4"))+
  labs(x = "Latitude (meters)", y = "Difference in mean habitat suitability (%)")+
  theme(axis.title.y = element_text(angle=90, vjust=2.125, face="bold", size=18),
        axis.text.y = element_text(size=14)) +
  theme(panel.grid = element_line(color = "darkgrey",
                                  size = 0.3,
                                  linetype = 2))+
  scale_x_continuous(n.breaks = 20)+
  theme(axis.text.x = element_text(size=14),axis.title.x = element_text(angle=0, vjust = 0, face="bold", size=18))+
  
  theme(panel.border = element_rect(size=0.7, colour = "black", fill=NA, linetype=1))+
  theme(legend.position="none")+
  ggtitle("P. fraudulenta, 2050_RCP-6.0")+
  coord_flip()

#210060
x.matrix210060 = matrix(data = xvals_dif210060, ncol = 1, byrow = TRUE)
y.matrix = matrix(data = yvals_present, ncol = 1, byrow = TRUE)
df_210060 <- cbind(x.matrix210060,y.matrix)
class(try)
df_210060 <- as.data.frame(df_210060)
library(dplyr)
df_210060 <- df_210060 %>% 
  mutate(C = if_else(V1 > 0, "Gain", "Loss"))
df_210060$C <- as.factor(df_210060$C)

library(ggplot2)
plot210060 <- ggplot(data=df_210060, aes(x=V2, y=V1, fill = C)) +
  geom_bar(stat="identity")+
  geom_hline(yintercept = 0, linetype="dashed", 
             color = "black", size=1.2)+
  ylim(-0.022,0.022)+
  scale_fill_manual(values=c("red4",
                             "royalblue4"))+
  labs(x = "Latitude (meters)", y = "Difference in mean habitat suitability (%)")+
  theme(axis.title.y = element_text(angle=90, vjust=2.125, face="bold", size=18),
        axis.text.y = element_text(size=14)) +
  theme(panel.grid = element_line(color = "darkgrey",
                                  size = 0.3,
                                  linetype = 2))+
  scale_x_continuous(n.breaks = 20)+
  theme(axis.text.x = element_text(size=14),axis.title.x = element_text(angle=0, vjust = 0, face="bold", size=18))+
  
  theme(panel.border = element_rect(size=0.7, colour = "black", fill=NA, linetype=1))+
  theme(legend.position="none")+
  ggtitle("P. fraudulenta, 2100_RCP-6.0")+
  coord_flip()

tiff(filename= 'fraudulenta_habdif205060.tiff',width = 12*ppi, height = 8*ppi, res = ppi)
plot205060
dev.off()
tiff(filename= 'fraudulenta_habdif210060.tiff',width = 12*ppi, height = 8*ppi, res = ppi)
plot210060
dev.off()


#205085
x.matrix205085 = matrix(data = xvals_dif205085, ncol = 1, byrow = TRUE)
y.matrix = matrix(data = yvals_present, ncol = 1, byrow = TRUE)
df_205085 <- cbind(x.matrix205085,y.matrix)
class(try)
df_205085 <- as.data.frame(df_205085)
library(dplyr)
df_205085 <- df_205085 %>% 
  mutate(C = if_else(V1 > 0, "Gain", "Loss"))
df_205085$C <- as.factor(df_205085$C)

library(ggplot2)
plot205085 <- ggplot(data=df_205085, aes(x=V2, y=V1, fill = C)) +
  geom_bar(stat="identity")+
  geom_hline(yintercept = 0, linetype="dashed", 
             color = "black", size=1.2)+
  ylim(-0.022,0.022)+
  scale_fill_manual(values=c("red4",
                             "royalblue4"))+
  labs(x = "Latitude (meters)", y = "Difference in mean habitat suitability (%)")+
  theme(axis.title.y = element_text(angle=90, vjust=2.125, face="bold", size=18),
        axis.text.y = element_text(size=14)) +
  theme(panel.grid = element_line(color = "darkgrey",
                                  size = 0.3,
                                  linetype = 2))+
  scale_x_continuous(n.breaks = 20)+
  theme(axis.text.x = element_text(size=14),axis.title.x = element_text(angle=0, vjust = 0, face="bold", size=18))+
  
  theme(panel.border = element_rect(size=0.7, colour = "black", fill=NA, linetype=1))+
  theme(legend.position="none")+
  ggtitle("P. fraudulenta, 2050_RCP-8.5")+
  coord_flip()

#210085
x.matrix210085 = matrix(data = xvals_dif210085, ncol = 1, byrow = TRUE)
y.matrix = matrix(data = yvals_present, ncol = 1, byrow = TRUE)
df_210085 <- cbind(x.matrix210085,y.matrix)
class(try)
df_210085 <- as.data.frame(df_210085)
library(dplyr)
df_210085 <- df_210085 %>% 
  mutate(C = if_else(V1 > 0, "Gain", "Loss"))
df_210085$C <- as.factor(df_210085$C)

library(ggplot2)
plot210085 <- ggplot(data=df_210085, aes(x=V2, y=V1, fill = C)) +
  geom_bar(stat="identity")+
  geom_hline(yintercept = 0, linetype="dashed", 
             color = "black", size=1.2)+
  ylim(-0.022,0.022)+
  scale_fill_manual(values=c("red4",
                             "royalblue4"))+
  labs(x = "Latitude (meters)", y = "Difference in mean habitat suitability (%)")+
  theme(axis.title.y = element_text(angle=90, vjust=2.125, face="bold", size=18),
        axis.text.y = element_text(size=14)) +
  theme(panel.grid = element_line(color = "darkgrey",
                                  size = 0.3,
                                  linetype = 2))+
  scale_x_continuous(n.breaks = 20)+
  theme(axis.text.x = element_text(size=14),axis.title.x = element_text(angle=0, vjust = 0, face="bold", size=18))+
  
  theme(panel.border = element_rect(size=0.7, colour = "black", fill=NA, linetype=1))+
  theme(legend.position="none")+
  ggtitle("P. fraudulenta, 2100_RCP-8.5")+
  coord_flip()

tiff(filename= 'fraudulenta_habdif205085.tiff',width = 12*ppi, height = 8*ppi, res = ppi)
plot205085
dev.off()
tiff(filename= 'fraudulenta_habdif210085.tiff',width = 12*ppi, height = 8*ppi, res = ppi)
plot210085
dev.off()

library(gridExtra)
tiff(filename= 'fraudulenta_habdif.tiff',width = 16*ppi, height = 20*ppi, res = ppi)
grid.arrange(plot205026,plot210026,
             plot205045,plot210045,
             plot205060,plot210060,
             plot205085,plot210085, nrow = 4)
dev.off()


#CENTROID PLOTS
#«Finally, additionalStats() calculates the centroid of the species range 
#(the mean latitude and longitude of the occupied pixels) for each time and scenario.»

centroid <- read.csv("C:\\Users\\franc\\OneDrive\\Área de Trabalho\\Academic\\SDM_ASP_paper\\Post_analysis\\centroid_analysis.csv")

str(centroid)

centroid$Species <- as.factor(centroid$Species)
centroid$RCP <- as.factor(centroid$RCP)
centroid$Time <- as.factor(centroid$Time)

#688834.3095 - present CentroidY fraudulenta
#1697659.914 - present CentroidY seriata
#1237476.996 - present CentroidY australis

View(centroid)
centroidfraudulenta <- subset(centroid, centroid$Species == "fraudulenta")
centroidfraudulenta <- subset(centroidfraudulenta, centroidfraudulenta$RCP !="present")
str(centroidfraudulenta)

library(ggplot2)
centroidfraudulentaplot<- ggplot(centroidfraudulenta, aes(x=RCP,CentroidY, fill = Time)) +
  geom_point(size=9, shape=21)+
  geom_hline(yintercept = 600000.1, linetype="dashed", 
             color = "royalblue4", size=1.5)+
labs(x = "", y = "")+
  theme(axis.title.y = element_text(angle=90, vjust=2.125, face="bold", size=18),
        axis.text.y = element_text(size=14)) +
  theme(panel.grid = element_line(color = "darkgrey",
                                  size = 0.3,
                                  linetype = 2))+
  theme(axis.text.x = element_text(size=14),axis.title.x = element_text(angle=0, vjust = 0, face="bold", size=18))+
  
  theme(panel.border = element_rect(size=0.7, colour = "black", fill=NA, linetype=1))+
  theme(legend.position="right")+
  ggtitle("P. fraudulenta")

centroidaustralis <- subset(centroid, centroid$Species == "australis")
centroidaustralis <- subset(centroidaustralis, centroidaustralis$RCP !="present")
str(centroidaustralis)

library(ggplot2)
centroidaustralisplot<- ggplot(centroidaustralis, aes(x=RCP,CentroidY, fill = Time)) +
  geom_point(size=9, shape=21)+
  geom_hline(yintercept = 525442.2, linetype="dashed", 
             color = "royalblue4", size=1.5)+
  labs(x = "", y = "Centroid Y Latitude (meters)")+
  theme(axis.title.y = element_text(angle=90, vjust=2.125, face="bold", size=18),
        axis.text.y = element_text(size=14)) +
  theme(panel.grid = element_line(color = "darkgrey",
                                  size = 0.3,
                                  linetype = 2))+
  theme(axis.text.x = element_text(size=14),axis.title.x = element_text(angle=0, vjust = 0, face="bold", size=18))+
  
  theme(panel.border = element_rect(size=0.7, colour = "black", fill=NA, linetype=1))+
  theme(legend.position="right")+
  ggtitle("P. australis")

centroidseriata <- subset(centroid, centroid$Species == "seriata")
centroidseriata <- subset(centroidseriata, centroidseriata$RCP !="present")
str(centroidseriata)

library(ggplot2)
centroidseriataplot<- ggplot(centroidseriata, aes(x=RCP,CentroidY, fill = Time)) +
  geom_point(size=9, shape=21)+
  geom_hline(yintercept = 1180409.5, linetype="dashed", 
             color = "royalblue4", size=1.5)+
  labs(x = "RCP", y = "")+
  theme(axis.title.y = element_text(angle=90, vjust=2.125, face="bold", size=18),
        axis.text.y = element_text(size=14)) +
  theme(panel.grid = element_line(color = "darkgrey",
                                  size = 0.3,
                                  linetype = 2))+
  theme(axis.text.x = element_text(size=14),axis.title.x = element_text(angle=0, vjust = 0, face="bold", size=18))+
  
  theme(panel.border = element_rect(size=0.7, colour = "black", fill=NA, linetype=1))+
  theme(legend.position="right")+
  ggtitle("P. seriata")

ppi = 360
library(ggpubr)
tiff(filename= 'centroids.tiff',width = 12*ppi, height = 8*ppi, res = ppi)
ggarrange(centroidaustralisplot,centroidseriataplot,centroidfraudulentaplot, nrow = 1,common.legend = TRUE, legend="top")
dev.off()









