change name to weightRaster

NEWS.md

@@ -1,6 +1,9 @@
 # roads 1.1.2
 * Change dynamic least cost paths (DLCP) to iterative least cost paths (ILCP) throughout
 * Change default `roadMethod` to `"ilcp"` in `projectRoads` 
+* Add ability to use a custom `weightFunction` and add a `weightFunction` `slopePenaltyFn` that determines the grade between two cells
+* Change argument name from `cost` to `weightRaster` since it no longer represents a cost surface
+and can now be inputs to the `weightFunction`. Also change `roadsInCost` to `roadsInWeight`.
 
 # roads 1.1.1
 * Fix an issue where updates to `terra` were causing roads to break

R/buildSimList.R

@@ -21,27 +21,27 @@
 #' or terra objects.
 #' 
 #' @param roads roads input
-#' @param cost cost input
+#' @param weightRaster weightRaster input
 #' @param roadMethod method of road projection
 #' @param landings landings input
-#' @param roadsInCost Whether the roads have already been burned into cost
+#' @param roadsInWeight Whether the roads have already been burned into the `weightRaster`.
 #' @param sim A sim list to update rather than building new.
 #' @noRd
 
-buildSimList <- function(roads, cost, roadMethod, landings, roadsInCost, 
+buildSimList <- function(roads, weightRaster, roadMethod, landings, roadsInWeight, 
                          sim = NULL){
-  if(!is(cost, "SpatRaster")){
-    if(is(cost, "RasterLayer")){
-      cost <- terra::rast(cost)
+  if(!is(weightRaster, "SpatRaster")){
+    if(is(weightRaster, "RasterLayer")){
+      weightRaster <- terra::rast(weightRaster)
     } else {
-      stop("cost must be provided as a SpatRaster or RasterLayer", call. = FALSE)
+      stop("weightRaster must be provided as a SpatRaster or RasterLayer", call. = FALSE)
     }
   } 
-  if(terra::minmax(cost)[1] == 0 || 0 %in% terra::unique(cost)[[1]]){
-    message("0s detected in cost raster, these will be considered as existing roads")
-  } else if(roadsInCost){
-    warning("No 0s detected in cost raster. If existing roads have not been ",
-            "included in the cost raster set roadsInCost = FALSE to have them ",
+  if(terra::minmax(weightRaster)[1] == 0 || 0 %in% terra::unique(weightRaster)[[1]]){
+    message("0s detected in weightRaster raster, these will be considered as existing roads")
+  } else if(roadsInWeight){
+    warning("No 0s detected in weightRaster. If existing roads have not been ",
+            "included in the weightRaster set roadsInWeight = FALSE to have them ",
             "burned in", call. = FALSE)
   }
 
@@ -50,10 +50,10 @@ buildSimList <- function(roads, cost, roadMethod, landings, roadsInCost,
   }
 
   # Burn roads into raster if not already for raster roads before converting 
-  if(!roadsInCost && is(roads, "SpatRaster")){
-    message("Burning in roads to cost raster from road raster")
-    cost <- cost * (roads == 0)
-    roadsInCost <- TRUE
+  if(!roadsInWeight && is(roads, "SpatRaster")){
+    message("Burning in roads to weightRaster raster from road raster")
+    weightRaster <- weightRaster * (roads == 0)
+    roadsInWeight <- TRUE
   }
 
   if(!(is(roads, "sf") || is(roads, "sfc"))){
@@ -121,9 +121,9 @@ buildSimList <- function(roads, cost, roadMethod, landings, roadsInCost,
         ) %>%
         sf::st_cast("POINT") %>% 
         sf::st_set_agr("constant") %>% 
-        sf::st_set_crs(sf::st_crs(cost))
+        sf::st_set_crs(sf::st_crs(weightRaster))
 
-      message("CRS of landings supplied as a matrix is assumed to match the cost")
+      message("CRS of landings supplied as a matrix is assumed to match the weightRaster")
     }else {
       stop("landings must be either SpatRaster, RasterLayer, sf object, SpatialPoints*, ",
            "or SpatialPolygons*",
@@ -144,49 +144,49 @@ buildSimList <- function(roads, cost, roadMethod, landings, roadsInCost,
 
   # check crs error if different
   if(!all(sf::st_crs(roads) == sf::st_crs(landings),
-          sf::st_crs(roads) == sf::st_crs(cost))){
-    stop("the crs of roads, landings and cost must match", call. = FALSE)
+          sf::st_crs(roads) == sf::st_crs(weightRaster))){
+    stop("the crs of roads, landings and weightRaster must match", call. = FALSE)
   }
 
-  # burn in roads to have 0 cost 
-  if(!roadsInCost){
-    message("Burning in roads to cost raster from sf")
+  # burn in roads to have 0 weight 
+  if(!roadsInWeight){
+    message("Burning in roads to weightRaster from sf")
 
-    cost <- burnRoadsInCost(roads, cost)
+    weightRaster <- burnRoadsInWeight(roads, weightRaster)
 
   }
 
-  # crop landings and roads to bbox of cost raster
+  # crop landings and roads to bbox of weightRaster raster
   nrland <- nrow(landings)
   nrroads <- nrow(roads)
 
   if(nrroads == 0){
     stop("nrow(roads) is 0. Please supply at least one existing road", call. = FALSE)
   }
 
-  ext <- sf::st_bbox(cost) %>% as.numeric() %>% 
+  ext <- sf::st_bbox(weightRaster) %>% as.numeric() %>% 
     `names<-`(c("xmin", "ymin", "xmax", "ymax"))
   landings <- sf::st_crop(sf::st_set_agr(landings, "constant"), ext) %>% sf::st_set_agr("constant")
   roads <- sf::st_crop(sf::st_set_agr(roads, "constant"), ext) %>% sf::st_set_agr("constant")
 
   if(nrland != nrow(landings)){
-    stop(nrland - nrow(landings), " landings are outside the cost surface. ",
-         "The cost surface must cover the extent of landings", call. = FALSE)
+    stop(nrland - nrow(landings), " landings are outside the weightRaster. ",
+         "The weightRaster must cover the extent of landings", call. = FALSE)
   }
 
   if(nrroads != nrow(roads)){
-    stop(nrroads - nrow(roads), " roads are outside the cost surface. ",
-         "The cost surface must cover the extent of roads", call. = FALSE)
+    stop(nrroads - nrow(roads), " roads are outside the weightRaster. ",
+         "The weightRaster must cover the extent of roads", call. = FALSE)
   }
 
   if(is.null(sim)){
-    sim <- rlang::env(roads = roads, costSurface = cost, 
+    sim <- rlang::env(roads = roads, weightRaster = weightRaster, 
                 roadMethod = roadMethod, 
                 landings = landings)
   } else {
     sim$roads <- roads
     sim$landings <- landings
-    sim$costSurface <- cost
+    sim$weightRaster <- weightRaster
     sim$roadMethod <- roadMethod
   }
   return(sim)
@@ -195,37 +195,37 @@ buildSimList <- function(roads, cost, roadMethod, landings, roadsInCost,
 
 
 
-#' Burn roads in to cost
+#' Burn roads in to weightRaster
 #'
-#' Use sf roads object to convert cost to 0 where roads already exist. This is
+#' Use sf roads object to convert weightRaster to 0 where roads already exist. This is
 #' an internal function and does not contain many checks. Use at own risk.
 #'
 #' @param roads sf object with road lines
-#' @param cost SpatRaster with cost of road development
+#' @param weightRaster SpatRaster with weights
 #'
-#' @return SpatRaster of cost with 0 for roads.
+#' @return a SpatRaster with 0 for roads.
 #' 
 #' @noRd
-burnRoadsInCost <- function(roads, cost){
+burnRoadsInWeight <- function(roads, weightRaster){
   # The crs is checked above but stars requires that they be identical
   if(!is.na(sf::st_crs(roads))){
-    roads <- sf::st_transform(roads, sf::st_crs(cost))
+    roads <- sf::st_transform(roads, sf::st_crs(weightRaster))
   }
 
   if(any(grepl("POINT", sf::st_geometry_type(roads, by_geometry = TRUE))) &&
      any(grepl("LINESTRING", sf::st_geometry_type(roads, by_geometry = TRUE)))){
     geom_types <- c("POINT", "LINESTRING")
 
-    rasts <- lapply(geom_types, function(x, rds, cst){
+    rasts <- lapply(geom_types, function(x, rds, wt){
       geom_roads <- sf::st_collection_extract(rds, type = x)
-      geom_rast <- terra::rasterize(terra::vect(geom_roads), cst,
+      geom_rast <- terra::rasterize(terra::vect(geom_roads), wt,
                                     background = 0) > 0
-    }, rds = roads, cst = cost)
+    }, rds = roads, wt = weightRaster)
 
     roadsRast <- !(rasts[[1]]|rasts[[2]])
   } else {
-    roadsRast <- terra::rasterize(terra::vect(roads), cost, background = 0) == 0 
+    roadsRast <- terra::rasterize(terra::vect(roads), weightRaster, background = 0) == 0 
   }
 
-  cost <- cost * roadsRast
+  weightRaster <- weightRaster * roadsRast
 }

R/buildSnapRoads.R

@@ -38,11 +38,11 @@ buildSnapRoads <- function(sim, roadsOut){
   # add to existing roads
   sim$roads <- bind_rows(sim$roads, snap_roads_lines)
 
-  # burn into cost as 0
-  sim$cost <- burnRoadsInCost(sim$roads, sim$costSurface)
+  # burn into weightRaster as 0
+  sim$weightRaster <- burnRoadsInWeight(sim$roads, sim$weightRaster)
 
   if(roadsOut == "raster"){
-    sim$roads <- sim$cost == 0
+    sim$roads <- sim$weightRaster == 0
   }
 
   return(invisible(sim))

R/getClosestRoad.R

@@ -31,7 +31,7 @@ getClosestRoad <- function(sim, ordering = "closest"){
   # lines to points at raster centers
   if(sf::st_geometry_type(roads.pts, by_geometry = FALSE) == "GEOMETRY"){
     rd.lines <- sf::st_collection_extract(roads.pts, type = "LINESTRING")
-    rd.lines.pts <- terra::extract(sim$costSurface, terra::vect(rd.lines), xy = TRUE) %>% 
+    rd.lines.pts <- terra::extract(sim$weightRaster, terra::vect(rd.lines), xy = TRUE) %>% 
       sf::st_as_sf(coords = c("x", "y"), crs = sf::st_crs(rd.lines)) %>% 
       sf::st_geometry()
 
@@ -53,7 +53,7 @@ getClosestRoad <- function(sim, ordering = "closest"){
 
   # find landings that are within the space of one raster cell from the road
   touching_road <- which(distToRoad < 
-                           terra::res(sim$costSurface)[1])
+                           terra::res(sim$weightRaster)[1])
 
   if(length(touching_road) > 0){
     # make snap roads for these ones 

R/getGraph.R

@@ -23,12 +23,12 @@
 #' @noRd
 
 getGraph<- function(sim, neighbourhood,method="old",weightFunction = function(x1,x2,...) (x1+x2)/2,...){
-  #sim = list(costSurface=costRaster);neighbourhood="octagon"
-  #gdistance method takes more time and less memory. See testAltGraphFns in RoadPaper repo for details.  resolution=res(sim$costSurface)[1]
+  #sim = list(weightRaster=weightRaster);neighbourhood="octagon"
+  #gdistance method takes more time and less memory. See testAltGraphFns in RoadPaper repo for details.  resolution=res(sim$weightRaster)[1]
 
   inargs = list(...)
   if(!is.element("resolution",names(inargs))){
-    resolution = terra::res(sim$costSurface)[1]
+    resolution = terra::res(sim$weightRaster)[1]
   }else{
     resolution = inargs$resolution
   }
@@ -46,7 +46,7 @@ getGraph<- function(sim, neighbourhood,method="old",weightFunction = function(x1
                   rook=4,
                   octagon=8,
                   queen=8)
-    x = gdistance::transition(as(sim$costSurface, "Raster"), transitionFunction=function(x) 1/weightFunction(x[1],x[2],resolution=resolution,...), directions=dirs)
+    x = gdistance::transition(as(sim$weightRaster, "Raster"), transitionFunction=function(x) 1/weightFunction(x[1],x[2],resolution=resolution,...), directions=dirs)
 
     if(neighbourhood=="octagon"){
       #correct for diagonal distances and other aspects of geographic distance
@@ -61,9 +61,9 @@ getGraph<- function(sim, neighbourhood,method="old",weightFunction = function(x1
     # define global varaiables to avoid NOTEs on R CMD check
     weightV <- from <- to <- w1 <- w2 <- NULL
 
-    # prepare the cost surface raster #===========
-    # get cost as data.table from raster
-    weightV <- data.table(weight = terra::values(sim$costSurface, mat = FALSE))
+    # prepare the weightRaster raster #===========
+    # get weight as data.table from raster
+    weightV <- data.table(weight = terra::values(sim$weightRaster, mat = FALSE))
 
     # get the id for ther verticies which is used to merge with the edge list from
     # adj
@@ -75,10 +75,10 @@ getGraph<- function(sim, neighbourhood,method="old",weightFunction = function(x1
       stop("neighbourhood type not recognized")
     }
 
-    nc <- terra::ncol(sim$costSurface)
-    ncel <- terra::ncell(sim$costSurface) %>% as.integer()
+    nc <- terra::ncol(sim$weightRaster)
+    ncel <- terra::ncell(sim$weightRaster) %>% as.integer()
 
-    edges <- terra::adjacent(sim$costSurface, cells = 1:as.integer(ncel),
+    edges <- terra::adjacent(sim$weightRaster, cells = 1:as.integer(ncel),
                              directions = "rook", pairs = TRUE) %>%
       data.table::as.data.table()
 
@@ -87,20 +87,20 @@ getGraph<- function(sim, neighbourhood,method="old",weightFunction = function(x1
 
     edges <- unique(edges)
 
-    # merge in the weights from a cost surface
+    # merge in the weights from a weightRaster
     edges.weight <- merge(x = edges, y = weightV, by.x = "from", by.y = "id")
 
     # reformat
     data.table::setnames(edges.weight, c("from", "to", "w1"))
 
-    # merge in the weights to a cost surface
+    # merge in the weights to a weightRaster
     edges.weight <- data.table::setDT(merge(x = edges.weight, y = weightV,
                                           by.x = "to", by.y = "id"))
 
     # reformat
     data.table::setnames(edges.weight, c("from", "to", "w1", "w2"))
 
-    # take the average cost between the two pixels and remove w1 w2
+    # apply the weightFunction across the two pixels and remove w1 w2
     edges.weight[,`:=`(weight = weightFunction(w1,w2,resolution=resolution,...),
                      w1 = NULL,
                      w2 = NULL)]
@@ -113,7 +113,7 @@ getGraph<- function(sim, neighbourhood,method="old",weightFunction = function(x1
       } else {mW = 1}
       weightV[, weight := weight*mW]
 
-      edges <- terra::adjacent(sim$costSurface, cells = 1:as.integer(ncel),
+      edges <- terra::adjacent(sim$weightRaster, cells = 1:as.integer(ncel),
                                directions = "bishop", pairs = TRUE) %>%
         data.table::as.data.table()
 
@@ -122,20 +122,20 @@ getGraph<- function(sim, neighbourhood,method="old",weightFunction = function(x1
 
       edges <- unique(edges)
 
-      # merge in the weights from a cost surface
+      # merge in the weights from a weightRaster
       edges_bishop <- merge(x = edges, y = weightV, by.x = "from", by.y = "id")
 
       # reformat
       data.table::setnames(edges_bishop, c("from", "to", "w1"))
 
-      # merge in the weights to a cost surface
+      # merge in the weights to a weightRaster
       edges_bishop <- data.table::setDT(merge(x = edges_bishop, y = weightV,
                                               by.x = "to", by.y = "id"))
 
       # reformat
       data.table::setnames(edges_bishop, c("from", "to", "w1", "w2"))
 
-      # take the average cost between the two pixels and remove w1 w2
+      # apply weightFunction across the two pixels and remove w1 w2
       edges_bishop[,`:=`(weight = weightFunction(w1,w2,resolution=resolution,...),
                          w1 = NULL,
                          w2 = NULL)]
@@ -148,7 +148,7 @@ getGraph<- function(sim, neighbourhood,method="old",weightFunction = function(x1
 
     }
 
-    # get rid of NAs caused by barriers. Drop the w1 and w2 costs.
+    # get rid of NAs caused by barriers. Drop the w1 and w2 columns.
     edges.weight <- na.omit(edges.weight)
 
     # set the ids of the edge list. Faster than using as.integer(row.names())

R/lcpList.R

@@ -18,13 +18,13 @@ lcpList<- function(sim){
     return(invisible(sim))
   }
   ##Get a list of of cell indexs for to and from points
-  paths.matrix <- cbind(terra::cellFromXY(sim$costSurface,
+  paths.matrix <- cbind(terra::cellFromXY(sim$weightRaster,
                                          sf::st_coordinates(sim$landings)), 
-                      terra::cellFromXY(sim$costSurface, sim$roads.close.XY))
+                      terra::cellFromXY(sim$weightRaster, sim$roads.close.XY))
 
   sim$paths.list <- split(paths.matrix, 1:nrow(paths.matrix))
 
   rm(paths.matrix)
 
   return(invisible(sim))
-}
+}

R/mstList.R

@@ -25,9 +25,9 @@ mstList<- function(sim){
     return(invisible(sim))
   }
   # get cell indexs for new road start and end
-  mst.v <- rbind(terra::cellFromXY(sim$costSurface, 
+  mst.v <- rbind(terra::cellFromXY(sim$weightRaster, 
                                               sf::st_coordinates(sim$landings)),
-                           terra::cellFromXY(sim$costSurface,
+                           terra::cellFromXY(sim$weightRaster,
                                               sim$roads.close.XY))
   mst.v <- as.vector(mst.v)
   paths.matrix <- unique(mst.v) 
@@ -37,10 +37,10 @@ mstList<- function(sim){
     # get an adjaceny matrix given the cell numbers
     mst.adj <- igraph::distances(sim$g, paths.matrix, paths.matrix) 
 
-    # set the verticies names as the cell numbers in the costSurface
+    # set the verticies names as the cell numbers in the weightRaster
     rownames(mst.adj) <- paths.matrix
 
-    # set the verticies names as the cell numbers in the costSurface
+    # set the verticies names as the cell numbers in the weightRaster
     colnames(mst.adj) <- paths.matrix 
 
     # create a graph