R/funcs.R

#' Calculates the proportion of a species that could migrate from a cell to its
#' neighbours
#'
#' Calculates the proportion of a species that could migrate from a cell to its
#' neighbours given a raster, the maximum dispersal distance, and the number of
#' times the time-lapse will have
#' @param Raster a raster with the space where the species will be inhabiting
#' @param Distance the maximum dispersal distance of the species
#' @param Time the number of time-slices to be used
#' @return a dataframe with the cell id from, the cell id to, and beta, that is
#' the proportion of individuals that will go from cell from to cell to
#' @examples
#' data("r")
#' DistConect(Raster = r, Distance = 1000000, Time = 7)
#' @importFrom dplyr filter
#' @importFrom dplyr group_by
#' @importFrom dplyr summarize
#' @importFrom gdistance accCost
#' @importFrom gdistance geoCorrection
#' @importFrom gdistance transition
#' @importFrom magrittr "%>%"
#' @importFrom raster ncell
#' @importFrom raster xyFromCell
#' @author Derek Corcoran <[email protected]>
#' @author Javier Fajardo <[email protected] >
#' @export

DistConect<- function(Raster, Distance, Time = 7){
  #First we make a transition layer with the function transition from gdistance
  h16  <- transition(Raster, transitionFunction=function(x){1},16,symm=FALSE)
  #Then geocorrect for projection
  h16   <- geoCorrection(h16, scl=FALSE)
  #Since transition layers work with XY rather than IDs, get a matrix of XY coordinates
  B <- xyFromCell(Raster, cell = 1:ncell(Raster))
  #This nested loop is where the Bottle neck is
  #Start a list
  connections <- list()
  #For each pair of cells in B
  accCost2 <- function(x, fromCoords) {

    fromCells <- cellFromXY(x, fromCoords)
    tr <- transitionMatrix(x)
    tr <- rBind(tr, rep(0, nrow(tr)))
    tr <- cBind(tr, rep(0, nrow(tr)))
    startNode <- nrow(tr)
    adjP <- cbind(rep(startNode, times = length(fromCells)), fromCells)
    tr[adjP] <- Inf
    adjacencyGraph <- graph.adjacency(tr, mode = "directed", weighted = TRUE)
    E(adjacencyGraph)$weight <- 1/E(adjacencyGraph)$weight
    return(shortest.paths(adjacencyGraph, v = startNode, mode = "out")[-startNode])
  }
  for (i in 1:nrow(B)){
    #Create a temporal raster for each row with the distance from cell xy to all other cells
    temp <- accCost2(h16,B[i,])
    index <- which(temp < Distance)
    connections[[i]] <- cbind(i, index, temp[index])
  }
  #Get everything together as a large data frame
  connections <- do.call("rbind", connections)
  connections <- as.data.frame(connections)
  colnames(connections) <- c("from", "to", "dist")
  connections$Beta <- exp(-(connections$dist/max(connections$dist)))
  b <- connections %>% group_by(from) %>% summarize(TotalBeta = sum(Beta))
  connections <-merge(connections, b)
  connections$beta <-connections$Beta /connections$TotalBeta
  connections<- dplyr::filter(connections, beta > quantile(beta, 0.05))
  connections <- connections[,c(1,2,6)]
  n <- nrow(connections)
  connections <- do.call("rbind", replicate((Time), connections, simplify = FALSE))
  connections$Time <- rep(c(0:(Time-1)), each =n)
  return(connections)
}

#' Generates an AMPL dat file from a Stack
#'
#' Generates an AMPL dat file from a Stack in which each file is the projection
#' of a species distribution model into a time slice
#' @param Stack a raster with the space where the species will be inhabiting
#' @param maxalpha the maximum rate of change of the population in optimal
#' coditions, defaults in 10
#' @param maxbiomass the maximum initial biomass of the population in optimal
#' coditions, defaults in 2
#' @param maxcapacidad the maximum biomass of the population in optimal coditions,
#' defaults in 10
#' @param name the name of the .dat file that will be exported, defaults in
#' stack
#' @param Dist the maximum dispersal distance of the species modeled in the
#' stack
#' @return exports a .dat file to feed the AMPL model
#' @examples
#' data("univariate")
#' RasterToAmplDat(Stack = univariate, Threshold = 0.5)
#' @importFrom raster nlayers
#' @importFrom raster values
#' @importFrom tidyr spread
#' @importFrom tidyr unite_
#' @author Derek Corcoran <[email protected]>
#' @author Javier Fajardo <[email protected] >
#' @export

RasterToAmplDat <- function(Stack, maxalpha = 10, maxbiomass = 2, maxcapacidad = 10, name = "Stack", Dist = 1000000, Threshold, costlayer){

  TempStack <- Stack
  values(TempStack)[values(TempStack) < Threshold] = 0
  values(TempStack)[values(TempStack) >= Threshold] = 1

  TempRaster <- sum(TempStack)

  TempRaster[values(TempRaster) > 0] = 1
  TempRaster[values(TempRaster) == 0] = NA

  Stack <- Stack*TempRaster

  Alpha <- list()
  for (i in 1:nlayers(Stack)){
    temp <- data.frame(Alpha = values(Stack[[i]]), ID = 1:length(values(Stack[[i]])), Time = i)
    Alpha[[i]] <- temp[complete.cases(temp),]
  }

  Alpha <-  do.call(rbind, Alpha)
  Alpha$Alpha <- Alpha$Alpha/max(Alpha$Alpha)*maxalpha


  Biomasa <- data.frame(Biomasa = (values(Stack[[1]])/max(values(Stack[[1]]), na.rm = T))*maxbiomass, ID = 1:length(values(Stack[[1]])))
  Biomasa <-  Biomasa[complete.cases(Biomasa),]
  Biomasa$Biomasa <- round(Biomasa$Biomasa, 4)

  Capacidad <- list()
  for (i in 1:nlayers(Stack)){
    temp <- data.frame(Capacidad = values(Stack[[i]]), ID = 1:length(values(Stack[[i]])), Time = i)
    Capacidad[[i]] <- temp[complete.cases(temp),]
  }

  Capacidad <-  do.call(rbind, Capacidad)
  Capacidad$Capacidad <- Capacidad$Capacidad/max(Capacidad$Capacidad)*maxcapacidad

  Nodos <- merge(Alpha, Capacidad)

  Alphas <- spread(Alpha, key = Time, value = Alpha)
  Alphas$ID <- paste("[", Alphas$ID, ",*]", sep = "")
  Alphas$T1 <- 0
  Alphas$T2 <- 1
  Alphas$T3 <- 2
  Alphas$T4 <- 3
  Alphas$T5 <- 4
  Alphas$T6 <- 5
  Alphas$T7 <- 6
  Alphas$T8 <- 7
  Alphas <- Alphas[,c(1,10,2,11,3,12,4,13,5,14,6,15,7,16,8,17,9)]
  Alphas$line <- "\n"

  Biomasas <- Biomasa[,c(2,1)]
  Biomasas$line <- "\n"

  Capacidades <- spread(Capacidad, key = Time, value = Capacidad)
  Capacidades$ID <- paste("[", Capacidades$ID, ",*]", sep = "")
  Capacidades$T1 <- 0
  Capacidades$T2 <- 1
  Capacidades$T3 <- 2
  Capacidades$T4 <- 3
  Capacidades$T5 <- 4
  Capacidades$T6 <- 5
  Capacidades$T7 <- 6
  Capacidades$T8 <- 7
  Capacidades <- Capacidades[,c(1,10,2,11,3,12,4,13,5,14,6,15,7,16,8,17,9)]
  Capacidades$line <- "\n"

  Cost <- data.frame(ID = Biomasas$ID, Cost = values(costlayer)[Biomasas$ID], line = "\n")

  Beta <- DistConect(Stack[[1]], Distance = Dist, Time = nlayers(Stack))
  Beta <- dplyr::filter(Beta, from %in% unique(Biomasa$ID),
                        to %in% unique(Biomasa$ID))
  temp <-  split(Beta, Beta$Time)
  Betas <- do.call(cbind, lapply(1:length(temp), function(i){
    if (i == 1){
      setNames(data.frame(paste("[",temp[[i]][["from"]], ",", temp[[i]][["to"]], ",*","]", sep = ""), temp[[i]]["Time"], temp[[i]]["beta"]),
               c("V", paste("T", i, sep = ""), i-1))
    } else if (i == length(temp)){
      setNames(data.frame(temp[[i]]["Time"], temp[[i]]["beta"], rep("\n", NROW(temp[[i]]))),
               c(paste("T", i, sep = ""), i-1, "line"))
    } else {
      setNames(data.frame(temp[[i]]["Time"], temp[[i]]["beta"]),
               c(paste("T", i, sep = ""), i-1))
    }
  }))

  sink(paste0(name, ".dat"))
  cat(c("set V :=", unique(Alpha$ID), ";"))
  cat("\n")
  cat("\n")
  cat(c("set E :=", paste0("(",unique(unite_(Beta, col = "V", sep = ",", from = c("from", "to"))$V), ")"), ";"))
  cat("\n")
  cat("param T:= 7;")
  cat("\n")
  cat("param alpha :=")
  cat("\n")
  cat(do.call(paste, Alphas))
  cat(";")
  cat("\n")
  cat("param beta :=")
  cat("\n")
  cat(do.call(paste, Betas))
  cat(";")
  cat("\n")
  cat("param b0 :=")
  cat("\n")
  cat(do.call(paste, Biomasas))
  cat(";")
  cat("\n")
  cat("param u :=")
  cat("\n")
  cat(do.call(paste, Capacidades))
  cat(";")
  cat("\n")
  cat("param c :=")
  cat("\n")
  cat(do.call(paste, Cost))
  cat(";")
  cat("\n")
  sink()
  return(list(Nodos = Nodos, Biomasa = Biomasa, Alphas = Alphas, Alpha = Alpha, Cost = Cost))
}

#' Generates an AMPL dat file from a Stack
#'
#' Generates an AMPL dat file from a Stack in which each file is the projection
#' of a species distribution model into a time slice
#' @param Stack a raster with the space where the species will be inhabiting
#' @param Distance the maximum dispersal distance of the species modeled in the
#' stack
#' @param Time number of time slices present in the Stack
#' @param name the name of the .dat file that will be exported, defaults in
#' stack
#' @param Threshold minimum value in the model to allow the species to exist
#' @param costlayer raster with the costs of each cell
#' @return exports a .dat file to feed the AMPL model
#' @examples
#' data("univariate")
#' RtoQuadAmplDat(Stack = univariate, Distance = 1000000, Threshold = 0.5,
#' name = "TeSt")
#' @importFrom gdistance accCost
#' @importFrom gdistance geoCorrection
#' @importFrom gdistance transition
#' @importFrom raster ncell
#' @importFrom raster nlayers
#' @importFrom raster values
#' @importFrom raster xyFromCell
#' @importFrom tidyr unite_
#' @author Derek Corcoran <[email protected]>
#' @author Javier Fajardo <[email protected] >
#' @export


RtoQuadAmplDat <- function(Stack, Distance, Threshold, name){
  accCost2 <- function(x, fromCoords) {

    fromCells <- cellFromXY(x, fromCoords)
    tr <- transitionMatrix(x)
    tr <- rBind(tr, rep(0, nrow(tr)))
    tr <- cBind(tr, rep(0, nrow(tr)))
    startNode <- nrow(tr)
    adjP <- cbind(rep(startNode, times = length(fromCells)), fromCells)
    tr[adjP] <- Inf
    adjacencyGraph <- graph.adjacency(tr, mode = "directed", weighted = TRUE)
    E(adjacencyGraph)$weight <- 1/E(adjacencyGraph)$weight
    return(shortest.paths(adjacencyGraph, v = startNode, mode = "out")[-startNode])
  }

  values(Stack)[values(Stack) < Threshold] = 0
  values(Stack)[values(Stack) >= Threshold] = 1

  Suitability <- list()
  for (i in 1:nlayers(Stack)){
    temp <- data.frame(Suitability = values(Stack[[i]]), ID = 1:length(values(Stack[[i]])), Time = i-1)
    Suitability[[i]] <- temp[complete.cases(temp),]
  }

  Suitability <- do.call("rbind", Suitability)

  Raster <- sum(Stack)

  Raster[values(Raster) > 0] = 1
  Raster[values(Raster) == 0] = NA

  h16  <- transition(Raster, transitionFunction=function(x){1},16,symm=FALSE)

  h16   <- geoCorrection(h16, scl=FALSE)

  ID <-c(1:ncell(Raster))[!is.na(values(Raster))]

  B <- xyFromCell(Raster, cell = ID)

  connections <- list()
  #For each pair of cells in B
  for (i in 1:nrow(B)){
    #Create a temporal raster for each row with the distance from cell xy to all other cells
    temp <- accCost2(h16,B[i,])
    index <- which(temp < Distance)
    connections[[i]] <- cbind(ID[i], index, temp[index])
  }
  #Get everything together as a large data frame
  connections <- do.call("rbind", connections)
  connections <- as.data.frame(connections)
  colnames(connections) <- c("from", "to", "dist")
  temp <-  split(Suitability, Suitability$Time)
  Suitability <- do.call(cbind, lapply(1:length(temp), function(i){
    if (i == 1){
      setNames(data.frame(paste("[",temp[[i]][["ID"]], ",*","]", sep = ""), temp[[i]]["Time"], temp[[i]]["Suitability"]),
               c("V", paste("T", i, sep = ""), i-1))
    } else if (i == length(temp)){
      setNames(data.frame(temp[[i]]["Time"], temp[[i]]["Suitability"], rep("\n", NROW(temp[[i]]))),
               c(paste("T", i, sep = ""), i-1, "line"))
    } else {
      setNames(data.frame(temp[[i]]["Time"], temp[[i]]["Suitability"]),
               c(paste("T", i, sep = ""), i-1))
    }
  }))

  sink(paste0(name, ".dat"))
  cat(c("set V :=", unique(unique(connections$to), unique(connections$to)), ";"))
  cat("\n")
  cat("\n")
  cat(c("set E :=", paste0("(",unique(unite_(connections, col = "V", sep = ",", from = c("from", "to"))$V), ")"), ";"))
  cat("\n")
  cat("\n")
  cat(paste0("param T:= ", nlayers(Stack),";"))
  cat("\n")
  cat("\n")
  cat("param u :=")
  cat("\n")
  cat(do.call(paste, Suitability))
  cat(";")
  cat("\n")
  sink()
  return(list(connections = connections, Suitability = Suitability))
}
#' Generates an AMPL dat file from a list of Stacks
#'
#' Generates an AMPL dat file from a Stack in which each file is the projection
#' of a species distribution model into a time slice
#' @param Stacklist a list of Stacks with the space where the species will be
#' inhabiting
#' @param Dist the maximum dispersal distance of the species modeled in the
#' stack
#' @param name the name of the .dat file that will be exported
#' @param nchains the number of chains to go through
#' @param costlayer raster with the costs of each cell
#' @return exports a .dat file to feed the AMPL model
#' @examples
#' data("BinSpp")
#' MultiSppQuad(Stacklist = BinSpp, Dist = 1000000, name = "Two")
#' @importFrom dplyr filter
#' @importFrom gdistance accCost
#' @importFrom gdistance geoCorrection
#' @importFrom gdistance transition
#' @importFrom raster ncell
#' @importFrom raster nlayers
#' @importFrom raster values
#' @importFrom raster xyFromCell
#' @importFrom tidyr unite_
#' @author Derek Corcoran <[email protected]>
#' @author Javier Fajardo <[email protected] >
#' @export

MultiSppQuad <- function(Stacklist, Dist, name, nchains = 100, costlayer){

  Masklayer <- costlayer
  values(Masklayer) <- ifelse(is.na(values(Masklayer)), NA, 1)
  for (i in 1:length(Stacklist)){
    Stacklist[[i]] <- Stacklist[[i]] * Masklayer
  }
  accCost2 <- function(x, fromCoords) {

    fromCells <- cellFromXY(x, fromCoords)
    tr <- transitionMatrix(x)
    tr <- rBind(tr, rep(0, nrow(tr)))
    tr <- cBind(tr, rep(0, nrow(tr)))
    startNode <- nrow(tr)
    adjP <- cbind(rep(startNode, times = length(fromCells)), fromCells)
    tr[adjP] <- Inf
    adjacencyGraph <- graph.adjacency(tr, mode = "directed", weighted = TRUE)
    E(adjacencyGraph)$weight <- 1/E(adjacencyGraph)$weight
    return(shortest.paths(adjacencyGraph, v = startNode, mode = "out")[-startNode])
  }
  Suitabilities <- list()
  for(j in 1:length(Stacklist)){
    Suitability <- list()
    for (i in 1:nlayers(Stacklist[[j]])){
      temp <- data.frame(Suitability = values(Stacklist[[j]][[i]]), ID = 1:length(values(Stacklist[[j]][[i]])), Time = i-1)
      Suitability[[i]] <- temp[complete.cases(temp),]
    }
    Suitabilities[[j]]<- do.call("rbind", Suitability)
    Suitabilities[[j]]$Spp <- names(Stacklist)[j]
  }
  Suitability <- do.call("rbind", Suitabilities)
  s <- Suitability %>% group_by(ID) %>% summarise(SUMA = sum(Suitability)) %>% filter(SUMA > 0)
  Suitability <- Suitability[Suitability$ID %in% s$ID,]


  Spps <- unique(Suitability$Spp)

  Suitability <- Suitability[,c(4,1,2,3)]

  Suitabilities <- list()
  for (i in Spps){
    Suitabilities[[i]] <- dplyr::filter(Suitability, Spp == i)

    temp <-  split(Suitabilities[[i]], Suitabilities[[i]]$Time)
    Suitabilities[[i]] <- do.call(cbind, lapply(1:length(temp), function(i){
      if (i == 1){
        setNames(data.frame(paste("[",temp[[i]][["Spp"]],",",temp[[i]][["ID"]], ",*","]", sep = ""), temp[[i]]["Time"], temp[[i]]["Suitability"]),
                 c("V", paste("T", i, sep = ""), i-1))
      } else if (i == length(temp)){
        setNames(data.frame(temp[[i]]["Time"], temp[[i]]["Suitability"], rep("\n", NROW(temp[[i]]))),
                 c(paste("T", i, sep = ""), i-1, "line"))
      } else {
        setNames(data.frame(temp[[i]]["Time"], temp[[i]]["Suitability"]),
                 c(paste("T", i, sep = ""), i-1))
      }
    }))
  }

  Suitability <-do.call("rbind", Suitabilities)

  conns <- list()
  for(j in 1:length(Stacklist)){
    Raster <- sum(Stacklist[[j]])

    Raster[values(Raster) > 0] = 1
    Raster[values(Raster) == 0] = NA

    h16  <- transition(Raster, transitionFunction=function(x){1},16,symm=FALSE)

    h16   <- geoCorrection(h16, scl=FALSE)

    ID <-c(1:ncell(Raster))[!is.na(values(Raster))]

    B <- xyFromCell(Raster, cell = ID)

    connections <- list()
    #For each pair of cells in B
    for (i in 1:nrow(B)){
      #Create a temporal raster for each row with the distance from cell xy to all other cells
      temp <- accCost2(h16,B[i,])
      index <- which(temp < Dist)
      connections[[i]] <- cbind(ID[i], index, temp[index])
    }
    #Get everything together as a large data frame
    connections <- do.call("rbind", connections)
    connections <- as.data.frame(connections)
    colnames(connections) <- c("from", "to", "dist")
    connections$Sp <- names(Stacklist)[j]
    conns[[j]] <- connections
  }
  connections <- conns
  connections <- do.call("rbind", connections)

  Nchains <- data.frame(Spp = Spps, Nchains = nchains, Space = "\n")
  Cost <- values(costlayer)[unique(unique(connections$to), unique(connections$to))]


  sink(paste0(name, ".dat"))
  cat(c("set V :=", unique(unique(connections$to), unique(connections$to)), ";"))
  cat("\n")
  cat("\n")
  cat(c("set SP :=", names(Stacklist), ";"))
  cat("\n")
  cat("\n")
  cat(c("set c :=", Cost, ";"))
  cat("\n")
  cat("\n")
  cat(c("set E :=", paste0("(",unique(unite_(connections, col = "V", sep = ",", from = c("from", "to"))$V), ")"), ";"))
  cat("\n")
  cat("\n")
  cat(paste0("param T:= ", (nlayers(Stacklist[[1]])-1),";"))
  cat("\n")
  cat("\n")
  cat("param u :=")
  cat("\n")
  cat(do.call(paste, Suitability))
  cat(";")
  cat("\n")
  cat("param nchains := ")
  cat("\n")
  cat(do.call(paste, Nchains))
  cat(";")
  cat("\n")
  sink()
  return(list(connections = connections, Suitability = Suitability))
}
derek-corcoran-barrios/RtoAmpl documentation built on May 14, 2019, 10:33 a.m.