Project road locations based on existing roads, planned landings, and a weights raster that in conjunction with a weighting function determines the cost to build a road between two cells in the raster.
Usage
projectRoads(
landings = NULL,
weightRaster = NULL,
roads = NULL,
roadMethod = "ilcp",
plotRoads = FALSE,
mainTitle = "",
neighbourhood = "octagon",
weightFunction = function(x1, x2, ...) (x1 + x2)/2,
sim = NULL,
roadsOut = NULL,
roadsInWeight = TRUE,
ordering = "closest",
...
)
# S4 method for ANY,ANY,ANY,ANY,ANY,ANY,ANY,ANY,missing
projectRoads(
landings = NULL,
weightRaster = NULL,
roads = NULL,
roadMethod = "ilcp",
plotRoads = FALSE,
mainTitle = "",
neighbourhood = "octagon",
weightFunction = function(x1, x2, ...) (x1 + x2)/2,
sim = NULL,
roadsOut = NULL,
roadsInWeight = TRUE,
ordering = "closest",
...
)
# S4 method for ANY,ANY,ANY,ANY,ANY,ANY,ANY,ANY,list
projectRoads(
landings = NULL,
weightRaster = NULL,
roads = NULL,
roadMethod = "ilcp",
plotRoads = FALSE,
mainTitle = "",
neighbourhood = "octagon",
weightFunction = function(x1, x2, ...) (x1 + x2)/2,
sim = NULL,
roadsOut = NULL,
roadsInWeight = TRUE,
ordering = "closest",
...
)Arguments
- landings
sf polygons or points, 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.
- 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 = FALSEand they will be burned in.- roads
sf lines, SpatialLines*, RasterLayer, SpatRaster. Existing road network.
- roadMethod
Character. Options are "ilcp", "mst", "lcp", "snap".
- plotRoads
Boolean. Should the resulting road network be plotted. Default FALSE.
- mainTitle
Character. A title for the plot
- neighbourhood
Character. 'rook','queen', or 'octagon'. The cells that should be considered adjacent. 'octagon' option is a modified version of the queen's 8 cell neighbourhood in which diagonals weights are 2^0.5x higher than horizontal/vertical weights.
- weightFunction
function. Method for calculating the weight of an edge between two nodes from the value of the weights raster at each of those nodes (x1 and x2). Default is the mean. Functions should be symmetric, meaning that the value returned does not depend on the ordering of x1 and x2. All functions must include the arguments
x1,x2and....- sim
list. Returned from a previous iteration of
projectRoads. weightRaster, roads, androadMethodare ignored if asimlist is provided.- roadsOut
Character. Either "raster", "sf" or NULL. If "raster" roads are returned as a raster in the
simlist. If "sf" the roads are returned as an sf object which will contain lines if the roads input was sf lines but a geometry collection of lines and points if the roads input was a raster. The points in the geometry collection represent the existing roads while new roads are created as lines. If NULL (default) then the returned roads are sf if the input is sf or Spatial* and raster if the input was a raster.- roadsInWeight
Logical. The default is TRUE which means the
weightRasteris assumed to include existing roads as 0. If FALSE then the roads will be "burned in" to theweightRasterwith a weight of 0.- ordering
character. The order in which roads should be built to landings when
roadMethod = "ilcp". Options are "closest" (default) where landings closest to existing roads are accessed first, or "none" where landings are accessed in the order they are provided in.- ...
Optional additional arguments to weightFunction
Value
a list with components:
roads: the projected road network, including new and input roads.
weightRaster: the weights raster, updated to have 0 for new roads that were added.
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
weightRaster. This is updated based on the new roads so that vertices that were connected by new roads now have a weight of 0. This can be used to avoid recomputing the graph in a simulation with multiple time steps.
Details
Four different methods for projecting road networks have been implemented:
"snap": Connects each landing directly to the closest road without reference to the weights or other landings.
"lcp": Least Cost Path connects each landing to the closest point on the road by determining the least cost path based on the weights raster provided, it does not consider other landings.
"ilcp": Iterative Least Cost Path, same as "lcp" but it builds each path sequentially so that later roads will use earlier roads. The sequence of landings is determined by
orderingand is "closest" by default, the other option is "none" which will use the order that landings are supplied in."mst": Minimum Spanning Tree connects all landings to the road by determining the least cost path to the road or other landings based on the weight raster.
Examples
CLUSexample <- prepExData(CLUSexample)
doPlots <- interactive()
projectRoads(CLUSexample$landings, CLUSexample$cost, CLUSexample$roads,
"lcp", plotRoads = doPlots, mainTitle = "CLUSexample")
#> 0s detected in weightRaster raster, these will be considered as existing roads
#> $landings
#> Simple feature collection with 4 features and 0 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 0.5 ymin: 0.5 xmax: 4.5 ymax: 2.5
#> Geodetic CRS: WGS 84
#> geometry
#> 1 POINT (0.5 2.5)
#> 2 POINT (2.5 2.5)
#> 3 POINT (1.5 0.5)
#> 4 POINT (4.5 0.5)
#>
#> $roads
#> class : SpatRaster
#> dimensions : 5, 5, 1 (nrow, ncol, nlyr)
#> resolution : 1, 1 (x, y)
#> extent : 0, 5, 0, 5 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84
#> source(s) : memory
#> name : lyr.1
#> min value : FALSE
#> max value : TRUE
#>
#> $g
#> IGRAPH 8bee29e U-W- 25 72 --
#> + attr: weight (e/n)
#> + edges from 8bee29e:
#> [1] 1-- 2 1-- 6 1-- 7 2-- 3 2-- 6 2-- 7 2-- 8 3-- 4 3-- 7 3-- 8
#> [11] 3-- 9 4-- 5 4-- 8 4-- 9 4--10 5-- 9 5--10 6-- 7 6--11 6--12
#> [21] 7-- 8 7--11 7--12 7--13 8-- 9 8--12 8--13 8--14 9--10 9--13
#> [31] 9--14 9--15 10--14 10--15 11--12 11--16 11--17 12--13 12--16 12--17
#> [41] 12--18 13--14 13--17 13--18 13--19 14--15 14--18 14--19 14--20 15--19
#> [51] 15--20 16--17 16--21 16--22 17--18 17--21 17--22 17--23 18--19 18--22
#> [61] 18--23 18--24 19--20 19--23 19--24 19--25 20--24 20--25 21--22 22--23
#> [71] 23--24 24--25
#>
#> $weightRaster
#> class : SpatRaster
#> dimensions : 5, 5, 1 (nrow, ncol, nlyr)
#> resolution : 1, 1 (x, y)
#> extent : 0, 5, 0, 5 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84
#> source(s) : memory
#> name : lyr.1
#> min value : 0.00000
#> max value : 19.84622
#>
#> $roadMethod
#> [1] "lcp"
#>
# More realistic examples that take longer to run
# \donttest{
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")
#> 0s detected in weightRaster raster, these will be considered as existing roads
### 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")
#> 0s detected in weightRaster raster, these will be considered as existing roads
### 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")
#> 0s detected in weightRaster raster, these will be considered as existing roads
#> CRS of landings supplied as a matrix is assumed to match the weightRaster
## 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")
#> 0s detected in weightRaster raster, these will be considered as existing roads
# }