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 86b1757..2bc453f 100644 --- a/R/getLandingsFromTarget.R +++ b/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 diff --git a/R/mstList.R b/R/mstList.R index f4eb314..43791f3 100644 --- a/R/mstList.R +++ b/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) } diff --git a/R/projectRoads.R b/R/projectRoads.R index a31fbce..9b5fee8 100644 --- a/R/projectRoads.R +++ b/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 diff --git a/man/getLandingsFromTarget.Rd b/man/getLandingsFromTarget.Rd index 1311d94..0024562 100644 --- a/man/getLandingsFromTarget.Rd +++ b/man/getLandingsFromTarget.Rd @@ -7,10 +7,18 @@ getLandingsFromTarget(harvest, landingDens = NULL, sampleType = "centroid") } \arguments{ +<<<<<<< HEAD \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 and 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.} +>>>>>>> b96540ec721aa63598cf23506b7ccfc903e0b443 \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 e70f0e8..5ce663e 100644 --- a/man/projectRoads.Rd +++ b/man/projectRoads.Rd @@ -58,10 +58,16 @@ projectRoads( ) } \arguments{ +<<<<<<< HEAD \item{landings}{sf polygons or points, RasterLayer, SpatialPolygons*, SpatialPoints*, or matrix, containing features to be connected +======= +\item{landings}{sf polygons or points, SpatRaster, RasterLayer, SpatialPolygons*, +SpatialPoints*, matrix, containing features to be connected +>>>>>>> b96540ec721aa63598cf23506b7ccfc903e0b443 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. A \code{weightRaster} and \code{weightFunction} together determine the cost to build a road between two adjacent raster cells. diff --git a/tests/testthat/test-getGraph.R b/tests/testthat/test-getGraph.R index 6c65357..f284721 100644 --- a/tests/testthat/test-getGraph.R +++ b/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", { diff --git a/tests/testthat/test-getLandingsFromTarget.R b/tests/testthat/test-getLandingsFromTarget.R index 700efce..dc0648f 100644 --- a/tests/testthat/test-getLandingsFromTarget.R +++ b/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) + }) diff --git a/tests/testthat/test-projectRoads.R b/tests/testthat/test-projectRoads.R index 3ea2693..49b9ae8 100644 --- a/tests/testthat/test-projectRoads.R +++ b/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 diff --git a/vignettes/articles/getDistFromSource.Rmd b/vignettes/articles/getDistFromSource.Rmd index 7396f8d..5244c7e 100644 --- a/vignettes/articles/getDistFromSource.Rmd +++ b/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. diff --git a/vignettes/grade-penalty.Rmd b/vignettes/grade-penalty.Rmd index f76eb9e..b7fd919 100644 --- a/vignettes/grade-penalty.Rmd +++ b/vignettes/grade-penalty.Rmd @@ -17,19 +17,5 @@ knitr::opts_chunk$set( ``` ```{r setup} -library(roads) - -demoScen <- prepExData(demoScen) - -purrr::map(demoScen, "cost.rast") %>% terra::rast() %>% terra::plot() - -can_elev <- geodata::elevation_30s("CAN", "test") - -ext <- terra::draw() - -ex_elev <- terra::crop(can_elev, ext) - -plot(ex_elev) - -# Move the above to a data-raw file +# see data being prepared in dem_example.R ``` diff --git a/vignettes/roads-vignette.Rmd b/vignettes/roads-vignette.Rmd index 36a0275..e9cbb16 100644 --- a/vignettes/roads-vignette.Rmd +++ b/vignettes/roads-vignette.Rmd @@ -65,7 +65,7 @@ CLUSexample <- prepExData(CLUSexample) ### 1. Weights Raster and Weight Function -The cost of new roads development within the landscape is determined by applying the `weightFunction` to the `weightRaster`. The `weightRaster` must be provided as a single, numeric `SpatRaster` or `RasterLayer` object. The cell values in the raster for two adjacent cells are provided to the `weightFunction` to calculate the weight for the edge connecting the two cells in the graph which should represent the cost of building a road connecting the two cells. Two ways of calculating the cost are supported by the package with the option for users to develop their own. The simplest method is for the `weightRaster` values to represent the cost of building a road across a cell and the for the `weightFunction` to simply take the mean of the values in adjacent cells. The other option included in the package is to use the `gradePenaltyFn`. This function takes a modified DEM raster as input and determines costs by applying a base cost and a grade penalty which is multiplied by the % grade between adjacent cells. The DEM can be modified by adding barriers to the raster, either as negative values which represent high costs (e.g. a stream crossing) or NAs for inaccessible areas (e.g. a cliff or a lake). See `?gradePenaltyFn` for more details. +The cost of new roads development within the landscape is determined by applying the `weightFunction` to the `weightRaster`. The `weightRaster` must be provided as a single, numeric `SpatRaster` or `RasterLayer` object. The cell values in the raster for two adjacent cells and the distance between the cells are provided to the `weightFunction` to calculate the weight for the edge connecting the two cells in the graph which should represent the cost of building a road connecting the two cells. Two ways of calculating the cost are supported by the package with the option for users to develop their own. The simplest method is for the `weightRaster` values to represent the cost of building a road across a cell and the for the `weightFunction` to simply take the mean of the values in adjacent cells (`simpleCostFn`). The other option included in the package is to use the `gradePenaltyFn`. This function takes a modified DEM raster as input and determines costs by applying a base cost and a grade penalty which is multiplied by the % grade between adjacent cells. The DEM can be modified by adding barriers to the raster, either as negative values which represent high costs (e.g. a stream crossing) or NAs for inaccessible areas (e.g. a cliff or a lake). See `?gradePenaltyFn` for more details. The following points apply to all `weightRasters`: