resolving merge conflicts

Merge branch 'master' of https://github.com/LandSciTech/roads

# Conflicts:
#	R/getLandingsFromTarget.R
#	R/projectRoads.R
#	man/getLandingsFromTarget.Rd
#	man/projectRoads.Rd

NEWS.md

@@ -2,11 +2,14 @@
 
 * 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` `gradePenaltyFn` 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`.
+* Add ability to use a custom `weightFunction` and add a `weightFunction` `gradePenaltyFn` 
+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`.
 * returned roads are no longer unioned together. 
 * Deprecate `getDistFromSource` and use `terra::distance` instead. 
+* Fix bug in getLandingsFromTarget and change so that patches are used for 0,1 rasters and 
+ids are used otherwise using terra::as.polygons to make it faster. 
 
 # roads 1.1.1
 * Fix an issue where updates to `terra` were causing roads to break

R/buildSimList.R

@@ -81,34 +81,8 @@ buildSimList <- function(roads, weightRaster, roadMethod, landings, roadsInWeigh
         sf::st_set_agr("constant")
 
     } else if(is(landings, "Raster") || is(landings, "SpatRaster")){
-      if(is(landings, "Raster")){
-        landings <- terra::rast(landings)
-      }
-
-      if(terra::nlyr(landings) > 1){
-        stop("landings should be a single layer SpatRaster")
-      }
-
-      # check if landings are clumps of cells (ie polygons) or single cells
-      # (ie points) and if clumps take centroid
-      clumpedRast <- terra::patches(landings, allowGaps = FALSE, zeroAsNA = TRUE) 
-
-      clumps <- clumpedRast %>% 
-        terra::freq()
-
-      clumps <- clumps[,2] %>% max() > 1
-
-      if(clumps){
-        landings <- sf::st_as_sf(terra::as.polygons(clumpedRast, 
-                                                          dissolve = TRUE)) %>% 
-          sf::st_set_agr("constant")
-      } else {
-        landings <- terra::subst(landings, from = 0, to = NA)%>% 
-          terra::as.points() %>% 
-          sf::st_as_sf() %>% 
-          sf::st_set_agr("constant")
-      }
 
+      landings <- getLandingsFromTarget(landings)
 
     } else if(is(landings, "matrix")){
       xyind <- which(colnames(landings) %in% c("x", "X", "y", "Y"))

R/getLandingsFromTarget.R

@@ -29,10 +29,11 @@
 #' area is determined by the CRS. For projected CRS this should likely be a very
 #' small number i.e. < 0.001.
 #'
-#' @param harvest sf, SpatialPolygons or RasterLayer object with harvested
-#'   areas. If it is a RasterLayer with more than one unique value other than 0
-#'   each value will be run separately which will produce different results from
-#'   a 0/1 raster and will be much slower.
+#' @param harvest sf, SpatialPolygons, SpatRaster or RasterLayer object with harvested
+#'   areas. If it is a raster with values outside 0,1, values are assumed
+#'   to be harvest block IDs. If raster values are in 0,1 they are assumed to be
+#'   a binary raster and `terra::patches` is used to identify harvest
+#'   blocks.
 #' @param landingDens number of landings per unit area. This should be in the
 #'   same units as the CRS of the harvest. Note that 0.001 points per m2 is > 1000
 #'   points per km2 so this number is usually very small for projected CRS.
@@ -114,16 +115,34 @@ getLandingsFromTarget <- function(harvest,
         harvest <- terra::rast(harvest)
       }
 
-      # check if harvest are clumps of cells (ie polygons) or single cells
-      # (ie points) and if clumps take centroid
-      clumpedRast <- terra::patches(harvest, allowGaps = FALSE, zeroAsNA = TRUE)
+      if(terra::nlyr(harvest) > 1){
+        stop("landings should be a single layer SpatRaster")
+      }
+
+      #if harvest rast is binary use clumping else assume values are cutblock ids
+      lnds_mnmx <- terra::minmax(harvest)[,1]
+
+      if(all(lnds_mnmx %in% c(0,1))){
+        message("harvest raster values are all in 0,1. Using patches as landing areas")
+        # check if harvest are clumps of cells (ie polygons) or single cells
+        # (ie points) and if clumps take centroid
+        harvest <- terra::patches(harvest, allowGaps = TRUE, zeroAsNA = TRUE)
+
+      } else {
+        message("harvest raster values not in 0,1. Using values as landing ids")
+
+      }
+
+      clumps <- harvest %>% terra::freq()
 
-      clumps <- clumpedRast %>%
-        terra::freq()
-      clumps <- max(clumps[,3])  > 1
+      clumps <- clumps[,3] %>% max() > 1
 
       if(clumps){
-        landings <- getLandingsFromTargetRast(harvest, landingDens, sampleType)
+        # zeros should be treated as NA
+        harvest <- terra::subst(harvest, from = 0, to = NA)
+        harvest <- sf::st_as_sf(terra::as.polygons(harvest,
+                                                   aggregate = TRUE)) %>%
+          sf::st_set_agr("constant")
       } else {
         landings <- terra::subst(harvest, from = 0, to = NA)%>%
           terra::as.points() %>%
@@ -135,12 +154,11 @@ getLandingsFromTarget <- function(harvest,
                   " landingDens is ignored and cells > 0 converted to points",
                   call. = FALSE)
         }
-
+        return(landings)
       }
-      return(landings)
+
     }
   }
-
   if(sf::st_geometry_type(harvest, by_geometry = FALSE) %in%
      c("POLYGON", "MULTIPOLYGON")){
     if(sampleType == "centroid"){
@@ -182,7 +200,7 @@ getLandingsFromTarget <- function(harvest,
 
         landings <- landings1$lands
       }
-      sf::st_sf(landings)
+      return(sf::st_sf(landings))
     }
 
   }
@@ -241,7 +259,7 @@ getLandingsFromTargetRast<-function(inputPatches,
 
     nl <- ifelse(is.null(landingDens), 0, landingDens)
 
-    ip <- inputPatches == i
+    ip <- inputPatches == clumpVals[i]
 
     ip[ip == 0] <- NA
 

R/mstList.R

@@ -48,7 +48,7 @@ mstList<- function(sim, roadsConnected){
 
       mst.adj <- cbind(rbind(lnd_to_lnd,t(lnd_to_rd)), rbind(lnd_to_rd, 0))
     } else {
-      both.v <- c(lnd.v, rd.v)
+      both.v <- unique(c(lnd.v, rd.v))
       mst.adj <- igraph::distances(sim$g, both.v, both.v)
       dimnames(mst.adj) <- list(both.v, both.v)
     }

R/projectRoads.R

@@ -34,6 +34,7 @@
 #'   * "snap": Connects each landing to the closest (by Euclidean distance) road without,
 #'   reference to the weights or other landings.
 #'
+<<<<<<< HEAD
 #' @param landings sf polygons or points, RasterLayer, SpatialPolygons*,
 #'   SpatialPoints*, or matrix, containing features to be connected
 #'   to the road network. Matrix should contain columns x, y with coordinates,
@@ -53,6 +54,18 @@
 #'     interpreted as 100% grade).
 #' @param roads sf lines, SpatialLines*, RasterLayer, SpatRaster. The existing road network.
 #' @param roadMethod Character. Options are "ilcp", "mst", "lcp", and "snap".
+=======
+#' @param landings sf polygons or points, SpatRaster, RasterLayer, SpatialPolygons*,
+#'   SpatialPoints*, matrix, containing features to be connected
+#'   to the road network. Matrix should contain columns x, y with coordinates,
+#'   all other columns will be ignored. Polygon and raster inputs will be
+#'   processed by `getLandingsFromTarget` to get the centroid of harvest blocks.
+#' @param weightRaster SpatRaster or RasterLayer. weights Raster where existing
+#'   roads must be the only cells with a weight of 0. If existing roads do not
+#'   have 0 weight set `roadsInWeight = FALSE` and they will be burned in.
+#' @param roads sf lines, SpatialLines*, RasterLayer, SpatRaster. Existing road network.
+#' @param roadMethod Character. Options are "ilcp", "mst", "lcp", "snap".
+>>>>>>> b96540ec721aa63598cf23506b7ccfc903e0b443
 #' @param plotRoads Boolean. Should the resulting road network be plotted.
 #'   Default FALSE.
 #' @param mainTitle Character. A title for the plot
@@ -211,6 +224,18 @@ setMethod(
                         roadMethod = roadMethod,
                         landings = landings,
                         roadsInWeight = roadsInWeight)
+
+    # Check memeory requirements
+    rast_cells <- terra::ncell(weightRaster)
+    ram_needed_Mb <- 10^(-1.6+0.769*log10(rast_cells))
+    ram_avail_Mb <- terra::free_RAM()/1000
+
+    if(ram_needed_Mb > ram_avail_Mb){
+      warning("This road projection is expected to require ", ram_needed_Mb, 
+              " Mb of RAM but only ", ram_avail_Mb, "Mb is available. ",
+              "Consider closing other applications and/or running on a larger machine")
+    }
+
 # browser()
     sim$landingsIn <- sim$landings
 

tests/testthat/test-getGraph.R

@@ -35,8 +35,7 @@ test_that("getGraph works with gdistance method", {
   expect_length(igraph::edge_attr(gR_gD, "weight"), 40)
   expect_length(igraph::edge_attr(gQ_gD, "weight"), 72)
   expect_length(igraph::edge_attr(gO_gD, "weight"), 72)
-  # TODO: figure out why not as expected but not really using gdistance right now
-  # expect_true(all(igraph::edge_attr(gO_gD, "weight") >= igraph::edge_attr(gQ_gD, "weight")))
+  expect_true(all(igraph::edge_attr(gO_gD, "weight") >= igraph::edge_attr(gQ_gD, "weight")))
 })
 
 test_that("getGraph works with gradePenaltyFun", {

tests/testthat/test-getLandingsFromTarget.R

@@ -101,22 +101,29 @@ test_that("raster with clumps input works no ID",{
 })
 
 test_that("raster with clumps input works with ID",{
-  rast <- demoScen[[1]]$landings.poly %>% terra::vect() %>%
-    terra::rasterize(terra::rast(demoScen[[1]]$cost.rast), field = "ID") %>% 
+  rst <- demoScen[[1]]$landings.poly %>% terra::vect() %>%
+    terra::rasterize(demoScen[[1]]$cost.rast, field = "ID") %>% 
     terra::`crs<-`(value = "EPSG:5070")
 
-  # make sure that a single celled havest block will work with clumps
-  rast[10,10] <- 6
+  # make sure that a single celled harvest block will work with clumps
+  rst[10,10] <- 20
 
-  # Show effect of ID
-  rast[78:88, 4:5] <- 7
+  # Show effect of ID and check for ID not sequential
+  rst[78:88, 4:5] <- 30
+
+  rst[is.na(rst)] <- 0
 
-  outRastCent <- getLandingsFromTarget(rast)
-  outRastRand <- getLandingsFromTarget(rast, landingDens = 0.1, 
+  outRastCent <- getLandingsFromTarget(rst)
+  outRastRand <- getLandingsFromTarget(rst, landingDens = 0.1,
                                        sampleType = "random")
-  outRastReg <- getLandingsFromTarget(rast, landingDens = 0.1, 
+  outRastReg <- getLandingsFromTarget(rst, landingDens = 0.1,
                                       sampleType = "regular")
 
+  land_vals <- terra::extract(rst, terra::vect(outRastCent), ID = FALSE) %>% pull(ID)
+
+  # all unique raster values represented in landings
+  expect_length(setdiff(land_vals, terra::unique(rst) %>% pull(ID)), 0)
+
   expect_type(outRastCent, "list")
 
   if(interactive()){
@@ -129,5 +136,11 @@ test_that("raster with clumps input works with ID",{
     terra::plot(rast)
     plot(outRastReg, col = "red", add = T)
   }
+
+  # compare to supplying raster to projectRoads
+  prRastCent <- projectRoads(rst, demoScen[[1]]$cost.rast, demoScen[[1]]$road.line)
+
+  expect_equal(prRastCent$landings, outRastCent)
+
 })
 

tests/testthat/test-projectRoads.R

@@ -148,6 +148,21 @@ test_that("duplicate roads are not created", {
 
 })
 
+test_that("landings on road or multiple landings in same cell work", {
+  CLUSexample <- prepExData(CLUSexample)
+
+  CLUSexample$landings <- bind_rows(CLUSexample$landings, 
+                                    list(sf::st_point(c(1.5, 4.2)), sf::st_point(c(1.5, 0.6))) %>% 
+                                      sf::st_as_sfc() %>%
+                                      sf::st_as_sf(crs = sf::st_crs(CLUSexample$landings)) %>%
+                                      rename(geometry = x))
+  expect_type(
+    projectRoads(CLUSexample$landings, CLUSexample$cost, CLUSexample$roads, 
+                 roadMethod = "mst"),
+    "list")
+
+
+})
 
 if(FALSE){
   # checking memory allocations

vignettes/articles/getDistFromSource.Rmd

@@ -39,7 +39,9 @@ rds_rast <- terra::rasterize(terra::vect(rds),
 terr_dist <- terra::distance(rds_rast, target = 0)
 
 # Or you can use the line directly
-terr_dist_vect <- terra::distance(rds)
+terr_dist_vect <- terra::distance(terra::vect(rds))
+
+terra::plot(terr_dist)
 
 ```
 
@@ -56,7 +58,8 @@ The third method "pfocal2" uses a global moving window to calculate the distance
 Below we will explore different settings and how they affect the speed and memory usage as well as the resulting distance map. The code below creates a source raster showing a curving road and applies the `getDistToSource` function under different parameters.
 
 ```{r run_dist, cache=TRUE}
-#make example roads from scratch
+if(requireNamespace("pfocal", quietly = TRUE)){
+  #make example roads from scratch
 rds <- data.frame(x = 1:1000, y = cos(1:1000/100)*500) %>% 
   sf::st_as_sf(coords = c("x", "y")) %>% 
   sf::st_union() %>% 
@@ -95,11 +98,14 @@ doFun <- function(x){
 
 # iterate over all parameters
 purrr::walk(1:nrow(param_grid), doFun)
+}
+
+
 ```
 
 ```{r make_plots, fig.height= 10, fig.width = 8}
 
-
+if(requireNamespace("pfocal", quietly = TRUE)){
 # all maps 
 all_maps <- purrr::map(1:nrow(param_grid), 
                        ~tmap::qtm(param_grid$result[[.x]], raster.style = "cont",
@@ -110,11 +116,13 @@ all_maps <- purrr::map(1:nrow(param_grid),
                                   layout.title.snap.to.legend = FALSE))
 
 tmap::tmap_arrange(all_maps, ncol = 3)
+}
 ```
 
 The maps above show the affect of the different parameters with `maxDist` changing the area covered, `kwidth` changing the smoothness of the gradations for the "terra" and "pfocal" methods and causing grid artefacts for the "pfocal2" method.
 
 ```{r make_plots2, fig.height= 7, fig.width = 10}
+if(requireNamespace("pfocal", quietly = TRUE)){
 get_hist <- function(rast){
   f <- terra::hist(rast, plot = FALSE, breaks = 30)
   dat <- data.frame(counts= f$counts,breaks = f$mids)
@@ -128,13 +136,14 @@ purrr::map_dfr(1:nrow(param_grid),
   geom_bar(stat = "identity",fill='blue',alpha = 0.8)+
   xlab("Distance")+ ylab("Frequency")+
   facet_grid(method +maxDist ~ kwidth, labeller = label_both)
-
+}
 ```
 
 
 Looking at histograms of the raster values shows the larger steps between distance values as `kwidth` increases for the "terra" and "pfocal" methods. For the "pfocal2" method the frequency of distance values stays similar with higher values of `kwidth` but the maps above show that some spatial accuracy is lost due to the aggregation and disaggregation processes.
 
 ```{r make_plots3, fig.height= 7, fig.width = 10}
+if(requireNamespace("pfocal", quietly = TRUE)){
 param_grid %>% 
   ggplot( aes(x = method, y = time))+
   geom_point()+
@@ -144,6 +153,7 @@ param_grid %>%
   ggplot( aes(x = method, y = mem_alloc))+
   geom_point()+
   facet_grid(kwidth~maxDist, labeller = label_both)
+}
 ```
 
 The last two graphs show that the "pfocal2" method is fastest in all the cases shown, however for very high `maxDist` values this may not be the case. The "pfocal" method uses the most memory of all the methods but note that this is the total memory allocated and does not account for deallocated memory, meaning not all the memory allocated is used at once.