diff --git a/NEWS.md b/NEWS.md index eec7c07..8964110 100644 --- a/NEWS.md +++ b/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 diff --git a/R/buildSimList.R b/R/buildSimList.R index 5521be2..9b45b34 100644 --- a/R/buildSimList.R +++ b/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")) diff --git a/R/getLandingsFromTarget.R b/R/getLandingsFromTarget.R index fba0d14..20bfe63 100644 --- a/R/getLandingsFromTarget.R +++ b/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 diff --git a/R/projectRoads.R b/R/projectRoads.R index ffeb4f8..c787f4d 100644 --- a/R/projectRoads.R +++ b/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. diff --git a/man/getLandingsFromTarget.Rd b/man/getLandingsFromTarget.Rd index 3d9dfce..066f960 100644 --- a/man/getLandingsFromTarget.Rd +++ b/man/getLandingsFromTarget.Rd @@ -7,10 +7,11 @@ getLandingsFromTarget(harvest, landingDens = NULL, sampleType = "centroid") } \arguments{ -\item{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.} +\item{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 \code{terra::patches} is used to identify harvest +blocks.} \item{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 diff --git a/man/projectRoads.Rd b/man/projectRoads.Rd index d2c6250..9f59244 100644 --- a/man/projectRoads.Rd +++ b/man/projectRoads.Rd @@ -58,10 +58,11 @@ projectRoads( ) } \arguments{ -\item{landings}{sf polygons or points, RasterLayer, SpatialPolygons*, +\item{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 \code{getLandingsFromTarget} to get the centroid of harvest blocks.} \item{weightRaster}{SpatRaster or RasterLayer. weights Raster where existing roads must be the only cells with a weight of 0. If existing roads do not diff --git a/tests/testthat/test-getLandingsFromTarget.R b/tests/testthat/test-getLandingsFromTarget.R index 700efce..ad02baa 100644 --- a/tests/testthat/test-getLandingsFromTarget.R +++ b/tests/testthat/test-getLandingsFromTarget.R @@ -105,11 +105,11 @@ 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, @@ -117,6 +117,11 @@ test_that("raster with clumps input works with ID",{ 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) + })