inst/doc/ReconstructRoadHistory.R

In roads: Road Network Projection

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(terra)
library(roads)
library(sf)
library(dplyr)

## colours for displaying cost raster 
if(requireNamespace("viridis", quietly = TRUE)){
  # Use colour blind friendly palette if available
  rastColours <- c('grey50', viridis::viridis(30))
} else {
  rastColours <- c('grey50', terrain.colors(30))
}

## using demo scenario 1
demoScen <- prepExData(demoScen)
scen <- demoScen[[1]]



## ----fig.width=6,fig.height=5-------------------------------------------------
# make "real" roads by using ilcp method
# use all landings and build roads to closest first
land.pnts <- scen$landings.points[scen$landings.points$set %in% c(1:4),]

realRoads <- projectRoads(land.pnts, scen$cost.rast, scen$road.line,
                          roadMethod = "ilcp")

plot(scen$cost.rast, col = rastColours, main = 'Scenario')
plot(realRoads$roads, add = TRUE)
plot(land.pnts, add = TRUE, pch = 21, cex = 2, bg = 'white')
text(st_coordinates(land.pnts), labels = land.pnts$set, cex = 0.6, adj = c(0.5, 0.3),
     xpd = TRUE)

## ----fig.width=6,fig.height=5-------------------------------------------------

# Use sequence of landings to assign sequence of roads built

# burn in realRoads to have cost of 0 so that projected roads will follow them.
roadsRast <- terra::rasterize(terra::vect(realRoads$roads), scen$cost.rast, 
                              background = 0) == 0

# doesn't work if weight is 0 because areas with 0 weight are assumed to already be
# roads so no new ones will be built
scen$cost.rast <- scen$cost.rast * (roadsRast + 0.00001) 

# set pre-existing roads to year 0
scen$road.line$Year <- 0

# initialize sim list with first landings set 
multiTime_sim <- list(projectRoads(land.pnts[land.pnts$set == 1,], scen$cost.rast, 
                                               scen$road.line))

multiTime_sim[[1]]$roads <- multiTime_sim[[1]]$roads %>% 
  mutate(Year = ifelse(is.na(Year), 1, Year))

# iterate over landings sets using the sim list from the previous run as input
for (i in 2:4) {
  newSim <- projectRoads(sim =  multiTime_sim[[i-1]], landings = land.pnts[land.pnts$set == i,])
  
  newSim$roads <- newSim$roads %>% 
    mutate(Year = ifelse(is.na(Year), i, Year))
  
  multiTime_sim <- c(
    multiTime_sim,
    list(newSim)
  ) 
}

plot(multiTime_sim[[4]]$roads["Year"], lwd = 2)