From 34ec980eb7083673dc970cc99a436c047b686742 Mon Sep 17 00:00:00 2001 From: Josie Hughes Date: Wed, 13 Mar 2024 10:38:46 -0400 Subject: [PATCH] experimenting with gdistance method for getting the graph --- R/getGraph.R | 189 ++++++++++++++++++++++++++++----------------------- 1 file changed, 103 insertions(+), 86 deletions(-) diff --git a/R/getGraph.R b/R/getGraph.R index ce1751c..1506b52 100644 --- a/R/getGraph.R +++ b/R/getGraph.R @@ -22,115 +22,132 @@ #' @param neighbourhood neighbourhood type #' @noRd -getGraph<- function(sim, neighbourhood){ - # define global varaiables to avoid NOTEs on R CMD check - w1 <- w2 <- NULL - - sim$paths.v <- NULL - # prepare the cost surface raster #=========== - # get cost as data.table from raster - weight <- data.table(weight = terra::values(sim$costSurface, mat = FALSE)) - - # get the id for ther verticies which is used to merge with the edge list from - # adj - weight[, id := seq_len(.N)] - - # get the adjacency using SpaDES function adj #======= - # rooks case - if(!is.element(neighbourhood, c("rook", "octagon", "queen"))) { - stop("neighbourhood type not recognized") - } - - nc <- terra::ncol(sim$costSurface) - ncel <- terra::ncell(sim$costSurface) %>% as.integer() +getGraph<- function(sim, neighbourhood,method="gdistance"){ + #gdistance method takes more time and less memory. See testAltGraphFns in RoadPaper repo for details. + if(method=="gdistance"){ + if(!is.element(neighbourhood, c("rook", "octagon"))) { + stop("neighbourhood type not recognized") + } + + dirs = switch(neighbourhood, + rook=4, + octagon=8) + x = gdistance::transition(raster::raster(sim$costSurface), transitionFunction=function(x) 1/mean(x), directions=dirs) + + y = gdistance::transitionMatrix(x) + sim$g <- igraph::graph_from_adjacency_matrix(y, mode="undirected", weighted=TRUE) + igraph::E(sim$g)$weight <- 1/igraph::E(sim$g)$weight + return(invisible(sim)) + }else{ + # define global varaiables to avoid NOTEs on R CMD check + w1 <- w2 <- NULL + + sim$paths.v <- NULL + # prepare the cost surface raster #=========== + # get cost as data.table from raster + weight <- data.table(weight = terra::values(sim$costSurface, mat = FALSE)) + + # get the id for ther verticies which is used to merge with the edge list from + # adj + weight[, id := seq_len(.N)] + + # get the adjacency using SpaDES function adj #======= + # rooks case + if(!is.element(neighbourhood, c("rook", "octagon", "queen"))) { + stop("neighbourhood type not recognized") + } + + nc <- terra::ncol(sim$costSurface) + ncel <- terra::ncell(sim$costSurface) %>% as.integer() - edges <- terra::adjacent(sim$costSurface, cells = 1:as.integer(ncel), - directions = "rook", pairs = TRUE) %>% - data.table::as.data.table() + edges <- terra::adjacent(sim$costSurface, cells = 1:as.integer(ncel), + directions = "rook", pairs = TRUE) %>% + data.table::as.data.table() - # switch to and from where to > from so can remove duplicates - edges[edges$from < edges$to, ] <- edges[edges$from < edges$to, c('to','from')] + # switch to and from where to > from so can remove duplicates + edges[edges$from < edges$to, ] <- edges[edges$from < edges$to, c('to','from')] - edges <- unique(edges) + edges <- unique(edges) - # merge in the weights from a cost surface - edges_rook <- merge(x = edges, y = weight, by.x = "from", by.y = "id") + # merge in the weights from a cost surface + edges_rook <- merge(x = edges, y = weight, by.x = "from", by.y = "id") - # reformat - data.table::setnames(edges_rook, c("from", "to", "w1")) + # reformat + data.table::setnames(edges_rook, c("from", "to", "w1")) - # merge in the weights to a cost surface - edges_rook <- data.table::setDT(merge(x = edges_rook, y = weight, - by.x = "to", by.y = "id")) + # merge in the weights to a cost surface + edges_rook <- data.table::setDT(merge(x = edges_rook, y = weight, + by.x = "to", by.y = "id")) - # reformat - data.table::setnames(edges_rook, c("from", "to", "w1", "w2")) + # reformat + data.table::setnames(edges_rook, c("from", "to", "w1", "w2")) - # take the average cost between the two pixels and remove w1 w2 - edges_rook[,`:=`(weight = (w1 + w2)/ 2, - w1 = NULL, - w2 = NULL)] + # take the average cost between the two pixels and remove w1 w2 + edges_rook[,`:=`(weight = (w1 + w2)/ 2, + w1 = NULL, + w2 = NULL)] - if(neighbourhood == "rook"){ - edges.weight = edges_rook - } else { + if(neighbourhood == "rook"){ + edges.weight = edges_rook + } else { - # bishop's case - multiply weights by 2^0.5 - if(neighbourhood == "octagon"){ - mW = 2^0.5 - } else {mW = 1} - weight[, weight := weight*mW] + # bishop's case - multiply weights by 2^0.5 + if(neighbourhood == "octagon"){ + mW = 2^0.5 + } else {mW = 1} + weight[, weight := weight*mW] - edges <- terra::adjacent(sim$costSurface, cells = 1:as.integer(ncel), - directions = "bishop", pairs = TRUE) %>% - data.table::as.data.table() + edges <- terra::adjacent(sim$costSurface, cells = 1:as.integer(ncel), + directions = "bishop", pairs = TRUE) %>% + data.table::as.data.table() - # edges[from < to, c("from", "to") := .(to, from)] - edges[edges$from < edges$to, ] <- edges[edges$from < edges$to, c('to','from')] + # edges[from < to, c("from", "to") := .(to, from)] + edges[edges$from < edges$to, ] <- edges[edges$from < edges$to, c('to','from')] - edges <- unique(edges) + edges <- unique(edges) - # merge in the weights from a cost surface - edges_bishop <- merge(x = edges, y = weight, by.x = "from", by.y = "id") + # merge in the weights from a cost surface + edges_bishop <- merge(x = edges, y = weight, by.x = "from", by.y = "id") - # reformat - data.table::setnames(edges_bishop, c("from", "to", "w1")) + # reformat + data.table::setnames(edges_bishop, c("from", "to", "w1")) - # merge in the weights to a cost surface - edges_bishop <- data.table::setDT(merge(x = edges_bishop, y = weight, - by.x = "to", by.y = "id")) + # merge in the weights to a cost surface + edges_bishop <- data.table::setDT(merge(x = edges_bishop, y = weight, + by.x = "to", by.y = "id")) - # reformat - data.table::setnames(edges_bishop, c("from", "to", "w1", "w2")) + # reformat + data.table::setnames(edges_bishop, c("from", "to", "w1", "w2")) - # take the average cost between the two pixels and remove w1 w2 - edges_bishop[,`:=`(weight = (w1 + w2)/ 2, - w1 = NULL, - w2 = NULL)] + # take the average cost between the two pixels and remove w1 w2 + edges_bishop[,`:=`(weight = (w1 + w2)/ 2, + w1 = NULL, + w2 = NULL)] - # get the edges list #================== - edges.weight = rbind(edges_rook, edges_bishop) - } + # get the edges list #================== + edges.weight = rbind(edges_rook, edges_bishop) + } - # get rid of NAs caused by barriers. Drop the w1 and w2 costs. - edges.weight <- na.omit(edges.weight) + # get rid of NAs caused by barriers. Drop the w1 and w2 costs. + edges.weight <- na.omit(edges.weight) - # set the ids of the edge list. Faster than using as.integer(row.names()) - edges.weight[, id := seq_len(.N)] + # set the ids of the edge list. Faster than using as.integer(row.names()) + edges.weight[, id := seq_len(.N)] - # clean-up - rm(edges_rook, edges_bishop, edges, weight, mW, nc, ncel) + # clean-up + rm(edges_rook, edges_bishop, edges, weight, mW, nc, ncel) - # make the graph #================ - edge_mat <- as.matrix(edges.weight)[,1:2] - # browser() - # create the graph using to and from columns. Requires a matrix input - sim$g <- igraph::graph.edgelist(edge_mat, dir = FALSE) + # make the graph #================ + edge_mat <- as.matrix(edges.weight)[,1:2] + # browser() + # create the graph using to and from columns. Requires a matrix input + sim$g <- igraph::graph.edgelist(edge_mat, dir = FALSE) - # assign weights to the graph. Requires a matrix input - igraph::E(sim$g)$weight <- as.matrix(edges.weight)[,3] + # assign weights to the graph. Requires a matrix input + igraph::E(sim$g)$weight <- as.matrix(edges.weight)[,3] - #------clean up - rm(edges.weight) - return(invisible(sim)) + #------clean up + rm(edges.weight) + return(invisible(sim)) + } }