#Set the working directory 
setwd("C:/Users/Dani/Desktop/2007")
dir()
#Calculate genetic distance
gen<-read.table("gammast.txt", header=T, row.names=1) #Read table with gammast values 
gen[upper.tri(gen)] <- t(gen)[upper.tri(gen)] #transform to a square matrix
Dgen<-as.dist(gen) #Calculate genetic distance

library(raster) #Load raster package
#Read a raster file with the study area and soil uses
map <- raster("comb_999.grd")
#map[is.na(map)] <- 0.01 #if there are any NA values: transform NA values to low numeric values (0.01)

library(gdistance) #Load gdistance package
coor<-read.table("lonlat.txt", header=T) #Read table with the coordinates (longitude-latitude) of the patches
cP<-cbind(coor$x, coor$y) #x=longitude and y=latitude. Geogrphical points of the patches

#Calculate geographic distance
geoDist <- pointDistance(cP, longlat=TRUE)
geoDist <- as.dist(geoDist) 

#Aggregate the map; aggregation groups rectangular areas to create larger cells
map1<- aggregate(map,6)

#Create a transition matrix from the RasterLayer
tr <- transition(map1, transitionFunction=mean, directions=8) #We set 8 direction which means that all cells are connected to their 8 neighbours 

#Geographic correction
trC <- geoCorrection(tr, type="c", scl=TRUE, multpl=FALSE) #for LCP (Least-Cost Path)
trR <- geoCorrection(tr, type="r", scl=TRUE, multpl=FALSE) #for IBR (Isoaltion By Resistance)

#calculate both the least-cost distance and the resistance distance between points
cosDist <- costDistance(trC,cP) #Least-cost
resDist <- commuteDistance(trR, cP) #Resistence


#Correlation between genetic distance and geographic distance
cor(Dgen,geoDist)
#Correlation between genetic distance and least-cost distance
cor(Dgen,cosDist)
#Correlation between genetic distance and resistance distance
cor(Dgen,resDist)

#Mantel test for IBD (Isolation by Distance), LCP, IBR
library(vegan)
IBD <- mantel(geoDist,Dgen, permutations=119);IBD
LCP <- mantel(cosDist,Dgen, permutations=119);LCP
IBR <- mantel(resDist,Dgen, permutations=119);IBR
#permutations = 119 due to small sizes of distance matrices (5x5). Max. number of permutations allowed (to avoid duplications)


###Graphs###
#pdf("figure_3.pdf", width=9, height = 3.2)
png("figure_3.png", unit="in", 9,3.2, res=600)
par(family = "Times new roman")
layout(matrix(c(1,2,3), 1, 3, byrow = T))


#a#
#IBD
plot(geoDist, Dgen, pch=19, main="Isolation by Distance", las=1,
     ylab=expression(gamma[ST]), xlab="Geographic distance (m)")
abline(lm(Dgen~geoDist), col="grey",lty=2, lwd=2)
legend("topleft", 
       c(paste("r =", round(IBD$statistic,3)), paste("p =", round(IBD$signif,3))),
       ncol=1, bty="n", y.intersp=0.8, xjust=0)

mtext("A)",3,line=1.5,cex=1.0, xpd=T, adj=-.03)


#b#
#LCP
plot((cosDist), Dgen, pch=19, main="Least-cost Path", las=1,
     ylab=expression(gamma[ST]), xlab=("Least-cost distance"))
abline(lm(Dgen~cosDist), col="grey",lty=2, lwd=2)
legend("topleft", 
       c(paste("r =", round(LCP$statistic,3)), paste("p =", round(LCP$signif,3))),
       ncol=1, bty="n", y.intersp=0.8, xjust=0)

mtext("B)",3,line=1.5,cex=1.0, xpd=T, adj=-.03)

#c#
#IBR
plot((resDist), Dgen, pch=19, main="Isolation by Resistance", las=1,
     ylab=expression(gamma[ST]), xlab=("Resistance distance"))

abline(lm(Dgen~resDist), col="grey",lty=2, lwd=2)
legend("topleft", 
       c(paste("r =", round(IBR$statistic,3)), paste("p =", round(IBR$signif,3))),
       ncol=1, bty="n", y.intersp=0.8, xjust=0)

mtext("C)",3,line=1.5,cex=1.0, xpd=T, adj=-.03)


dev.off()


#Quickest paths#
#Calculate the speed of movement between adjacent cells
#LCP
adj1<- adjacent(map1, cells=1:ncell(map1), pairs=T, directions=8)
speed1<- trC
speed1[adj1]<- 6*exp(-3.5 * abs(trC[adj1]+0.05))
LCP<-geoCorrection(speed1)#a diagonal connection between cells takes longer to cross than a straight connection (van Etten 2017). 

#Coordinate (longitude, latitude) for each patch
FR1 <-c(-72.54561, -39.49427)
FR2 <-c(-72.54658, -39.49028)
FR3 <-c(-72.52903, -39.48622)
FR4 <-c(-72.52653, -39.48789)
FR5 <- c(-72.53903, -39.49092)


#Assign colors to each soil use
colors<-c("#8c510a",
       "#8c510a",
       "#c7eae5",
       "#01665e",
       "#f6e8c3",
       "#c7eae5",
       "#5ab4ac",
       "#d8b365")


#Path lines for LCP

#FR1 - FR2# It did not estimate
#FR1 - FR3#
FR1aFR3<-shortestPath(LCP, FR1, FR3, output="SpatialLines")
#FR1 - FR4#
FR1aFR4<-shortestPath(LCP, FR1, FR4, output="SpatialLines")
#FR1 - FR5# It did not estimate
#FR1aFR5<-shortestPath(LCP, FR1, FR5, output="SpatialLines")
#FR2 - FR3#
FR2aFR3<-shortestPath(LCP, FR2, FR3, output="SpatialLines")
#FR2 - FR4#
FR2aFR4<-shortestPath(LCP, FR2, FR4, output="SpatialLines")
#FR2 - FR5# It did not estimate
#FR2aFR5<-shortestPath(LCP, FR2, FR5, output="SpatialLines")
#FR3 - FR4# It did not estimate
#FR3aFR4<-shortestPath(LCP, FR3, FR4, output="SpatialLines")
#FR3 - FR5#
FR3aFR5<-shortestPath(LCP, FR3, FR5, output="SpatialLines")
#FR4 - FR5#
FR4aFR5<-shortestPath(LCP, FR4, FR5, output="SpatialLines")



### PLOT ###
#Assign colors to each soil use (same colors of Figure 1)

#pdf("Figure_4.pdf", height = 13, width = 7)
png("Figure_4_26july.png", unit="in", height = 5, width = 6, res= 600)

par(mar=c(4,4,1,0))

#a#
plot(map, xlab="", cex.axis=.8, 
     ylab="Latitude", col=colors, xlim=c(-72.58, -72.50), ylim=c(-39.52, -39.46)
     , las = 1, legend=FALSE, axes=T)


library(sf)

areas<-read_sf("C:/Users/Dani/Dropbox/Dani/Mapa_datos/areas.shp") #Load patch area shapefile
plot(areas, add=T,axis=F, col="gray", ylab="", xlab="") #Plot patch areas

mtext(side=1, text="Longitude", line=2)


#Black lines

lines(FR1aFR3, col=rgb(0,0,0,0.8), lwd=6)
text(FR1[1] - 10, FR1[2] - 10, "FR1")
text(FR3[1] + 10, FR3[2] + 10, "FR3")

lines(FR1aFR4, col=rgb(0,0,0,0.8), lwd=6)
text(FR1[1] - 10, FR1[2] - 10, "FR1")
text(FR4[1] + 10, FR4[2] + 10, "FR4")

lines(FR2aFR3, col=rgb(0,0,0,0.8), lwd=6)
text(FR2[1] - 10, FR2[2] - 10, "FR2")
text(FR3[1] + 10, FR3[2] + 10, "FR3")

lines(FR2aFR4, col=rgb(0,0,0,0.8), lwd=6)
text(FR2[1] - 10, FR2[2] - 10, "FR2")
text(FR4[1] + 10, FR4[2] + 10, "FR4")

lines(FR3aFR5, col=rgb(0,0,0,0.8), lwd=6)
text(FR3[1] - 10, FR3[2] - 10, "FR3")
text(FR5[1] + 10, FR5[2] + 10, "FR5")

lines(FR4aFR5, col=rgb(0,0,0,0.8), lwd=6)
text(FR4[1] - 10, FR4[2] - 10, "FR4")
text(FR5[1] + 10, FR5[2] + 10, "FR5")

dev.off()