fixes #33 for raster landings

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

@@ -30,10 +30,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 but 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.
@@ -108,24 +109,40 @@ getLandingsFromTarget <- function(harvest,
 
   if(!(is(harvest, "sf") || is(harvest, "sfc"))){
     if(is(harvest, "Spatial")){
-
+      
       harvest <- sf::st_as_sf(harvest) %>% sf::st_set_agr("constant")
-
+      
     } else if(is(harvest, "Raster") || is(harvest, "SpatRaster")){
       if(is(harvest, "Raster")){
         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)
 
-      clumps <- clumpedRast %>%
-        terra::freq() 
-      clumps <- max(clumps[,3])  > 1
+      if(terra::nlyr(landings) > 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 <- clumps[,3] %>% max() > 1
 
       if(clumps){
-        landings <- getLandingsFromTargetRast(harvest, landingDens, sampleType)
+        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() %>% 
@@ -137,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"){
@@ -181,16 +197,17 @@ getLandingsFromTarget <- function(harvest,
 
         landings <- c(landings1$lands, landings2Plus)
       } else {
-
+        
         landings <- landings1$lands
       }
-      sf::st_sf(landings)
+      return(sf::st_sf(landings))
     }
-
+    
   }
 }
 
 
+
 #' Select random landing locations within patches.
 #'
 #' @param inputPatches A SpatRaster. Harvested patches should have values

R/projectRoads.R

@@ -36,10 +36,11 @@
 #'   determining the least cost path to the road or other landings based on the
 #'   weight raster.
 #'
-#' @param landings sf polygons or points, RasterLayer, SpatialPolygons*,
+#' @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.
+#'   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.

tests/testthat/test-getLandingsFromTarget.R

@@ -105,18 +105,23 @@ test_that("raster with clumps input works with ID",{
     terra::rasterize(terra::rast(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
+  rast[10,10] <- 20
 
-  # Show effect of ID
-  rast[78:88, 4:5] <- 7
+  # Show effect of ID and check for ID not sequential
+  rast[78:88, 4:5] <- 30
 
   outRastCent <- getLandingsFromTarget(rast)
   outRastRand <- getLandingsFromTarget(rast, landingDens = 0.1, 
                                        sampleType = "random")
   outRastReg <- getLandingsFromTarget(rast, landingDens = 0.1, 
                                       sampleType = "regular")
 
+  land_vals <- terra::extract(rast, terra::vect(outRastCent), ID = FALSE) %>% pull(ID)
+
+  # all unique raster values represented in landings
+  expect_length(setdiff(land_vals, terra::unique(rast) %>% pull(ID)), 0)
+
   expect_type(outRastCent, "list")
 
   if(interactive()){
@@ -129,5 +134,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(rast, demoScen[[1]]$cost.rast, demoScen[[1]]$road.line)
+
+  expect_equal(prRastCent$landings, outRastCent)
+
 })