Title:	Road Network Projection
Version:	1.2.0
Date:	2024-06-26
URL:	https://github.com/LandSciTech/roads, https://landscitech.github.io/roads/
Description:	Iterative least cost path and minimum spanning tree methods for projecting forest road networks. The methods connect a set of target points to an existing road network using 'igraph' https://igraph.org to identify least cost routes. The cost of constructing a road segment between adjacent pixels is determined by a user supplied weight raster and a weight function; options include the average of adjacent weight raster values, and a function of the elevation differences between adjacent cells that penalizes steep grades. These road network projection methods are intended for integration into R workflows and modelling frameworks used for forecasting forest change, and can be applied over multiple time-steps without rebuilding a graph at each time-step.
License:	Apache License (≥ 2)
Encoding:	UTF-8
LazyData:	true
Imports:	dplyr, igraph (≥ 2.0.3), data.table, sf, units, rlang, methods, tidyselect, terra
RoxygenNote:	7.3.1
Suggests:	testthat (≥ 2.1.0), knitr, rmarkdown, viridis, tmap, bench, gdistance
VignetteBuilder:	knitr
Depends:	R (≥ 2.10)
Collate:	'CLUSexample.R' 'buildSimList.R' 'buildSnapRoads.R' 'demoScen.R' 'getClosestRoad.R' 'getDistFromSource.R' 'getGraph.R' 'weightFunctions.R' 'getLandingsFromTarget.R' 'lcpList.R' 'mstList.R' 'pathsToLines.R' 'projectRoads.R' 'rasterToLineSegments.R' 'shortestPaths.R' 'plotRoads.R' 'rasterizeLine.R' 'prepExData.R' 'roads-package.R' 'dem_example.R'
BugReports:	https://github.com/LandSciTech/roads/issues
NeedsCompilation:	no
Packaged:	2024-06-27 14:29:56 UTC; EndicottS
Author:	Sarah Endicott [aut, cre], Kyle Lochhead [aut], Josie Hughes [aut], Patrick Kirby [aut], Her Majesty the Queen in Right of Canada as represented by the Minister of the Environment [cph] (Copyright holder for included functions buildSimList, getLandingsFromTarget, pathsToLines, plotRoads, projectRoads, rasterizeLine, rasterToLineSegments), Province of British Columbia [cph] (Copyright holder for included functions getGraph, lcpList, mstList, shortestPaths, getClosestRoad, buildSnapRoads)
Maintainer:	Sarah Endicott <sarah.endicott@ec.gc.ca>
Repository:	CRAN
Date/Publication:	2024-06-27 14:50:02 UTC

roads: Road Network Projection

Description

Iterative least cost path and minimum spanning tree methods for projecting forest road networks. The methods connect a set of target points to an existing road network using igraph https://igraph.org to identify least cost routes. The cost of constructing a road segment between adjacent pixels is determined by a user supplied 'weightRaster' and a 'weightFunction'; options include the average of adjacent 'weightRaster' values, and a function of the elevation differences between adjacent cells that penalizes steep grades. These road network projection methods are intended for integration into R workflows and modelling frameworks used for forecasting forest change, and can be applied over multiple timesteps without rebuilding a graph at each timestep.

Author(s)

Maintainer: Sarah Endicott sarah.endicott@ec.gc.ca (ORCID)

Authors:

• Kyle Lochhead Kyle.Lochhead@gov.bc.ca

• Josie Hughes josie.hughes@ec.gc.ca

• Patrick Kirby

Other contributors:

• Her Majesty the Queen in Right of Canada as represented by the Minister of the Environment (Copyright holder for included functions buildSimList, getLandingsFromTarget, pathsToLines, plotRoads, projectRoads, rasterizeLine, rasterToLineSegments) [copyright holder]

• Province of British Columbia (Copyright holder for included functions getGraph, lcpList, mstList, shortestPaths, getClosestRoad, buildSnapRoads) [copyright holder]

See Also

Useful links:

• https://github.com/LandSciTech/roads

• https://landscitech.github.io/roads/

• Report bugs at https://github.com/LandSciTech/roads/issues

Data from the CLUS example

Description

From Kyle Lochhead and Tyler Muhly's CLUS road simulation example. SpatRaster files created with the terra package must be saved with terra::wrap() and need to be unwrapped before they are used. prepExData() does this.

Usage

data(CLUSexample)

Format

A named list with components:

• cost: a PackedSpatRaster representing road building cost.

• landings: an sf dataframe of points representing landing locations.

• roads: a PackedSpatRaster representing existing roads.

Examples

CLUSexample
prepExData(CLUSexample)

Grade penalty example data

Description

A list containing two rasters covering an area near Revelstoke, British Columbia, Canada. ex_elev is elevation data and ex_wat is the proportion of the cell that contains water. Both are subsets of data downloaded with the geodata package at 30 arc seconds resolution.SpatRaster files created with the terra package must be saved with terra::wrap() and need to be unwrapped before they are used. prepExData() does this.

Usage

data(dem_example)

Format

A named list with components:

• ex_elev: a PackedSpatRaster of elevation.

• ex_wat: a PackedSpatRaster of proportion water.

Details

Elevation data are primarily from Shuttle Radar Topography Mission (SRTM), specifically the hole-filled CGIAR-SRTM (90 m resolution) from https://srtm.csi.cgiar.org/.

Water data are derived from the ESA WorldCover data set at 0.3-seconds resolution. (License CC BY 4.0). See https://esa-worldcover.org/en for more information.

References

Zanaga, D., Van De Kerchove, R., De Keersmaecker, W., Souverijns, N., Brockmann, C., Quast, R., Wevers, J., Grosu, A., Paccini, A., Vergnaud, S., Cartus, O., Santoro, M., Fritz, S., Georgieva, I., Lesiv, M., Carter, S., Herold, M., Li, Linlin, Tsendbazar, N.E., Ramoino, F., Arino, O., 2021. ESA WorldCover 10 m 2020 v100. doi:10.5281/zenodo.5571936.

Examples

dem_example
prepExData(dem_example)

Demonstration set of 10 input scenarios

Description

A demonstration set of scenarios that can be used as input to projectRoads(). The data contains SpatRaster objects that must be wrapped to be stored. To unwrap them use prepExData()

Usage

data(demoScen)

Format

A list of sub-lists, with each sub-list representing an input scenario. The scenarios (sub-lists) each contain the following components:

• scen.number: An integer value representing the scenario number (generated scenarios are numbered incrementally from 1).

• road.rast: A logical PackedSpatRaster representing existing roads. TRUE is existing road. FALSE is not existing road.

• road.line: A sf object representing existing roads.

• cost.rast: A PackedSpatRaster representing the cost of developing new roads on a given cell.

• landings.points: A sf object representing landings sets and landing locations within each set. The data frame includes a field named 'set' which contains integer values representing the landings set that each point belongs to

• landings.stack: A PackedSpatRaster with multiple layers representing the landings and landings sets. Each logical layer represents one landings set. Values of TRUE are a landing in the given set. Values of FALSE are not.

• landings.poly: A sf object representing a single set of polygonal landings.

See Also

projectRoads

Examples

demoScen[[1]]
demoScen <- prepExData(demoScen)
demoScen[[1]]

Moving window approach to get distance from source

Description

This function is deprecated please use terra::distance(). Note that you need to set target = 0 to get distances from cells that are zero to cells that are non-zero.

Usage

getDistFromSource(src, maxDist, kwidth = 3, method = "terra", override = FALSE)

Arguments

Details

This function provides three different methods for calculating the distance of all points on a landscape from "source" locations. This is a computationally intensive process so the function arguments can be used to balance the tradeoffs between speed and accuracy. Note the pfocal versions are only available in the development version of the package.

The "terra" and "pfocal" methods use an iterative moving window approach and assign each cell a distance based on the number of times the moving window is repeated before it is included. This means that the moving window function is run many times but for a small window relative to the size of the raster. The maxDist argument determines the maximum distance calculated and affects the number of iterations of the moving window that are needed. kwidth is the radius of the moving window in number of cells, with larger values reducing the number of iterations needed but also reducing the granularity of the distances produced. The resulting distances will be in increments of kwidth

• the resolution of the raster. The total number of iterations is maxDist/ kwidth * resolution. The only difference in these methods is the underlying package used to do the moving window. The terra package has methods for handling large rasters by writing them to disk, while the pfocal package requires that the raster can be held in memory as a matrix.

The third method "pfocal2" uses a global moving window to calculate the distance to the source. This means that the moving window only needs to be applied once but the window size can be very large. In this case maxDist determines the total size of the window. kwidth can be used to reduce the number of cells included in the moving window by aggregating the source raster by a factor of kwidth. This will increase the speed of computation but will produce results with artefacts of the larger grid and which may be less accurate since the output raster is disaggregated using bilinear interpolation.

Value

A SpatRaster

Examples

CLUSexample <-  prepExData(CLUSexample)
# Deprecated
# getDistFromSource(CLUSexample$roads, 5, 2)

# Use terra::distance instead
terra::distance(CLUSexample$roads, target = 0)


 library(sf)
 library(terra)

#make example roads from scratch
rds <- data.frame(x = 1:1000/100, y = cos(1:1000/100)) %>%
  st_as_sf(coords = c("x", "y")) %>%
  st_union() %>%
  st_cast("LINESTRING")

rds_rast <- rasterize(vect(rds),
                      rast(nrows = 50, ncols = 50,
                           xmin = 0, xmax = 10,
                           ymin = -5, ymax = 5),
                      touches = TRUE)

terra::distance(rds_rast)

# or straight from the line
terra::distance(rds_rast, terra::vect(rds %>% st_set_crs(st_crs(rds_rast))))

Get landing target points within harvest blocks

Description

Generate landing points inside polygons representing harvested area. There are three different sampling types available: "centroid" (default) returns the centroid or a point inside the polygon if the centroid is not (see sf::st_point_on_surface()); "random" returns a random sample given landingDens see (sf::st_sample()); "regular" returns points on a regular grid with cell size sqrt(1/landingDens) that intersect the polygon, or centroid if no grid points fall within the polygon.

Usage

getLandingsFromTarget(harvest, landingDens = NULL, sampleType = "centroid")

Arguments

Details

Note that the landingDens is points per unit area where the unit of area is determined by the CRS. For projected CRS this should likely be a very small number i.e. < 0.001.

Value

an sf simple feature collection with an ID column and POINT geometry

Examples

doPlots <- interactive()
demoScen <- prepExData(demoScen)

polys <- demoScen[[1]]$landings.poly[1:2,]

# Get centroid
outCent <- getLandingsFromTarget(polys)

if(doPlots){
  plot(sf::st_geometry(polys))
  plot(outCent, col = "red", add = TRUE)
}

# Get random sample with density 0.1 points per unit area
outRand <- getLandingsFromTarget(polys, 0.1, sampleType = "random")

if(doPlots){
  plot(sf::st_geometry(polys))
  plot(outRand, col = "red", add = TRUE)
}

# Get regular sample with density 0.1 points per unit area
outReg <- getLandingsFromTarget(polys, 0.1, sampleType = "regular")

if(doPlots){
  plot(sf::st_geometry(polys))
  plot(outReg, col = "red", add = TRUE)
}

Grade penalty edge weight function

Description

Method for calculating the weight of an edge between two nodes from the value of the input raster at each of those nodes (x1 and x2), designed for a single DEM input. The method assumes an input weightRaster in which:

• NA indicates a road cannot be built

• Negative values are costs for crossing streams or other barriers that are crossable but expensive. Edges that link to barrier penalty (negative value) nodes are assigned the largest barrier penalty weight.

• Zero values are assumed to be existing roads.

• All other values are interpreted as elevation in the units of the raster map (so that a difference between two cells equal to the map resolution can be interpreted as 100% grade) This is a simplified version of the grade penalty approach taken by Anderson and Nelson (2004): The approach does not distinguish between adverse and favourable grades. Default construction cost values are from the BC interior appraisal manual. The approach ignores (unknown) grade penalties beside roads and barriers in order to avoid increased memory and computational burden associated with multiple input rasters.

Usage

gradePenaltyFn(
  x1,
  x2,
  hdistance,
  baseCost = 16178,
  limit = 20,
  penalty = 504,
  limitWeight = NA
)

Arguments

References

Anderson AE, Nelson J (2004) Projecting vector-based road networks with a shortest path algorithm. Canadian Journal of Forest Research 34:1444–1457. https://doi.org/10.1139/x04-030

Examples

gradePenaltyFn(0.5,0.51,1)
gradePenaltyFn(0.5,0.65,1)
# grade > 20% so NA
gradePenaltyFn(0.5,0.75,1)

Plot projected roads

Description

Plot the results of projectRoads()

Usage

plotRoads(sim, mainTitle, subTitle = paste0("Method: ", sim$roadMethod), ...)

Arguments

Value

Creates a plot using base graphics

Examples

CLUSexample <- prepExData(CLUSexample)
prRes <- projectRoads(CLUSexample$landings, CLUSexample$cost, CLUSexample$roads)
if(interactive()){
  plotRoads(prRes, "Title")
}

Prepare example data

Description

Prepare example data included in the package that contain wrapped terra objects. This applies terra::unwrap() recursively to the list provided so that all PackedSpatRasters are converted to SpatRasters.

Usage

prepExData(x)

Arguments

Value

The same list but with unwrapped SpatRasters

Examples

CLUSexample
prepExData(CLUSexample)

Project road network

Description

Project a road network that links target landings to existing roads. For all methods except "snap", a weightRaster and weightFunction together determine the cost to build a road between two adjacent raster cells.

Usage

projectRoads(
  landings = NULL,
  weightRaster = NULL,
  roads = NULL,
  roadMethod = "ilcp",
  plotRoads = FALSE,
  mainTitle = "",
  neighbourhood = "octagon",
  weightFunction = simpleCostFn,
  sim = NULL,
  roadsOut = NULL,
  roadsInWeight = TRUE,
  ordering = "closest",
  roadsConnected = FALSE,
  ...
)

## S4 method for signature 'ANY,ANY,ANY,ANY,ANY,ANY,ANY,ANY,missing'
projectRoads(
  landings = NULL,
  weightRaster = NULL,
  roads = NULL,
  roadMethod = "ilcp",
  plotRoads = FALSE,
  mainTitle = "",
  neighbourhood = "octagon",
  weightFunction = simpleCostFn,
  sim = NULL,
  roadsOut = NULL,
  roadsInWeight = TRUE,
  ordering = "closest",
  roadsConnected = FALSE,
  ...
)

## S4 method for signature 'ANY,ANY,ANY,ANY,ANY,ANY,ANY,ANY,list'
projectRoads(
  landings = NULL,
  weightRaster = NULL,
  roads = NULL,
  roadMethod = "ilcp",
  plotRoads = FALSE,
  mainTitle = "",
  neighbourhood = "octagon",
  weightFunction = simpleCostFn,
  sim = NULL,
  roadsOut = NULL,
  roadsInWeight = TRUE,
  ordering = "closest",
  roadsConnected = FALSE,
  ...
)

Arguments

Details

Four road network projection methods are:

• "lcp": The Least Cost Path method connects each landing to the closest road with a least cost path, without reference to other landings.

• "ilcp": The Iterative Least Cost Path method iteratively connects each landing to the closest road with a least cost path, so that the path to each successive landing can include roads constructed to access previous landings. The sequence of landings is determined by ordering and is "closest" by default. The alternative "none" option processes landings in the order supplied by the user.

• "mst": The Minimum Spanning Tree method connects landings to the existing road with a minimum spanning tree that does not require users to specify the order in which landings are processed.

• "snap": Connects each landing to the closest (by Euclidean distance) road without, reference to the weights or other landings.

Value

a list with components:

• roads: the projected road network, including new and input roads.

• weightRaster: the updated weightRaster in which new and old roads have value 0.

• roadMethod: the road simulation method used.

• landings: the landings used in the simulation.

• g: the graph that describes the cost of paths between each cell in the updated weightRaster. Edges between vertices connected by new roads have weight 0. g can be used to avoid the cost of rebuilding the graph in a simulation with multiple time steps.

Examples

CLUSexample <- prepExData(CLUSexample)
doPlots <- interactive()

projectRoads(CLUSexample$landings, CLUSexample$cost, CLUSexample$roads,
             "lcp", plotRoads = doPlots, mainTitle = "CLUSexample")


# More realistic examples that take longer to run


demoScen <- prepExData(demoScen)

### using:  scenario 1 / sf landings / iterative least-cost path ("ilcp")
# demo scenario 1
scen <- demoScen[[1]]

# landing set 1 of scenario 1:
land.pnts <- scen$landings.points[scen$landings.points$set==1,]

prRes <- projectRoads(land.pnts, scen$cost.rast, scen$road.line, "ilcp",
                         plotRoads = doPlots, mainTitle = "Scen 1: SPDF-LCP")

### using: scenario 1 / `SpatRaster` landings / minimum spanning tree ("mst")
# demo scenario 1
scen <- demoScen[[1]]

# the RasterLayer version of landing set 1 of scenario 1:
land.rLyr <- scen$landings.stack[[1]]

prRes <- projectRoads(land.rLyr, scen$cost.rast, scen$road.line, "mst",
                         plotRoads = doPlots, mainTitle = "Scen 1: Raster-MST")


### using: scenario 2 / matrix landings raster roads / snapping ("snap")
# demo scenario 2
scen <- demoScen[[2]]

# landing set 5 of scenario 2, as matrix:
land.mat  <- scen$landings.points[scen$landings.points$set==5,] |>
  sf::st_coordinates()

prRes <- projectRoads(land.mat, scen$cost.rast, scen$road.rast, "snap",
                      plotRoads = doPlots, mainTitle = "Scen 2: Matrix-Snap")

## using scenario 7 / Polygon landings raster / minimum spanning tree
# demo scenario 7
scen <- demoScen[[7]]
# rasterize polygonal landings of demo scenario 7:
land.polyR <- terra::rasterize(scen$landings.poly, scen$cost.rast)

prRes <- projectRoads(land.polyR, scen$cost.rast, scen$road.rast, "mst",
                         plotRoads = doPlots, mainTitle = "Scen 7: PolyRast-MST")

Convert raster to lines

Description

Converts rasters that represent lines into an sf object.

Usage

rasterToLineSegments(rast, method = "mst")

Arguments

Details

For method = "nearest" raster is first converted to points and then lines are drawn between the nearest points. If there are two different ways to connect the points that have the same distance both are kept which can cause doubled lines. USE WITH CAUTION. method = "mst" converts the raster to points, reclassifies the raster so roads are 0 and other cells are 1 and then uses projectRoads to connect all the points with a minimum spanning tree. This will always connect all raster cells and is slower but will not double lines as often. Neither method is likely to work for very large rasters

Value

an sf simple feature collection

Examples

CLUSexample <- prepExData(CLUSexample)
# works well for very simple roads
roadLine1 <- rasterToLineSegments(CLUSexample$roads)

# longer running more realistic examples

demoScen <- prepExData(demoScen)
# mst method works well in this case
roadLine2 <- rasterToLineSegments(demoScen[[1]]$road.rast)

# nearest method has doubled line where the two roads meet
roadLine3 <- rasterToLineSegments(demoScen[[1]]$road.rast, method = "nearest")

# The mst method can also produce odd results in some cases
roadLine4 <- rasterToLineSegments(demoScen[[4]]$road.rast)

Faster rasterize for lines

Description

Rasterize a line using stars because fasterize doesn't work on lines and rasterize is slow. Deprecated use terra::rasterize

Usage

rasterizeLine(sfLine, rast, value)

Arguments

Value

a RasterLayer where the value of cells that touch the line will be the row index of the line in the sf

Examples

CLUSexample <- prepExData(CLUSexample)
roadsLine <- sf::st_sf(geometry = sf::st_sfc(sf::st_linestring(
matrix(c(0.5, 4.5, 4.5, 4.51),
       ncol = 2, byrow = TRUE) 
)))

# Deprecated rasterizeLine(roadsLine, CLUSexample$cost, 0)   
# Use terra::rasterize
terra::rasterize(roadsLine, CLUSexample$cost, background = 0)

Simple cost edge weight function

Description

Calculates the weight of an edge between two nodes as the mean value of an input cost raster at each of those nodes (x1 and x2).

Usage

simpleCostFn(x1, x2, hdistance)

Arguments

Examples

simpleCostFn(0.5,0.7,1)