R/SecondTrailsOpening.R
In VincyaneBadouard/LoggingLab: An R Package to Simulate Forest Logging
Defines functions
secondtrailsopening
Documented in
secondtrailsopening
#' Secondary skidding trails layout
#'
#' @description Starting from the main skidding trails, draw secondary skidding
#' trails on the zones accessible to the machines, allowing to collect the
#' selected trees with the chosen machines ("winching" argument). The layout
#' is optimised to reduce the distance covered while respecting topographical
#' constraints and avoiding trees to protect.
#'
#' @param topography Digital terrain model (DTM) of the inventoried plot (LiDAR
#' or SRTM) (\code{\link{DTMParacou}}) (RasterLayer **with a crs in UTM**)
#'
#' @param plotmask Inventoried plot mask
#' (SpatialPolygonsDataFrame **with a crs in UTM**)
#'
#' @param maintrails Main trails defined at the entire harvestable area (sf
#' linestring **with a crs in UTM**)
#'
#' @param treeselectionoutputs A list with:
#' - your inventory with: "DistCriteria", "Slope", "SlopeCriteria",
#' "LoggingStatus", "Selected", "Up", "VolumeCumSum",
#' "ProbedHollowProba", "ProbedHollow" new columns
#' (see the outputs metadata in the vignette).
#'
#' - the objective volume
#'
#' - 6 sets of spatial points:
#' harvestable, selected, future and reserve, hollow and fuel wood trees
#'
#' @param harvestablepolygons Accessible zones of the inventoried plot
#' (default: \code{\link{harvestableareadefinition}}) (sfc_MULTIPOLYGON)
#'
#' @param machinepolygons Accessible zones for machines of the inventoried plot
#' (default: \code{\link{harvestableareadefinition}}) (sf polygons data.frame)
#'
#' @param maintrailsaccess Access point of maintrail for each PU (prospection
#' unit) (sf or sfc)
#'
#' @param plotslope Slopes (in radians) of the inventoried plot (with a
#' neighbourhood of 8 cells) (default:
#' \code{\link{HarvestableAreaOutputsCable}})
#' (RasterLayer **with a crs in UTM**)
#'
#'@param scenario Logging scenario among: "RIL1", "RIL2broken", "RIL2", "RIL3",
#' "RIL3fuel", "RIL3fuelhollow" or "manual"(character) (see the
#' vignette)
#'
#'@param winching
#' "0": no cable or grapple (trail to tree foot)
#' "1": only cable (default = 40m)
#' "2": grapple (default = 6m) + cable (grapple priority)
#' If grapple + cable (winching = "2") without fuel wood (fuel = "0")
#' recovery of the tree foot with grapple if possible (respected grapple
#' conditions) otherwise with cable with angle to the trail.
#' Avoidance of future/reserves if chosen.
#'
#'@param fuel Fuel wood exploitation: no exploitation = "0", exploitation of
#' damage and unused part of logged trees for fuel = "1", exploitation of
#' hollow trees, damage and and unused part of the log for fuel = "2".
#' If fuel wood harvesting is chosen, these trails are preliminary,
#' and the trails will be adjusted after tree felling for fuel wood recovery.
#' The mortality of the trails is therefore not recorded at this stage
#'
#'@param verbose Allow to provide messages from internal functions (boolean)
#'
#' @param advancedloggingparameters Other parameters of the logging simulator
#' \code{\link{loggingparameters}} (list)
#'
#' @return A list with :
#' - *inventory*: Updated inventory
#'
#' - *SmoothedTrails*: Smoothed secondary trails (MULTIPOLYGON with crs)
#'
#' - *TrailsDensity*: Secondary trails density (in m/ha)
#'
#' - *TrailsIdentity*: Information on sections of the trails (matrix) with:
#' - *LineID*:
#' - *LoggedTrees*: idTree of trees reached by the trails
#' - *TypeExpl*: type of winching
#'
#' - *MainTrailsAccess*: Random access point of main trail for each harvestable
#' zone (sfc_POINT with crs)
#'
#' - *RawSecondTrails*: Non-smoothed secondary trails (SpatialLines with crs)
#'
#' - *CostRasterAgg*: The cost raster (RasterLayer with crs)
#'
#' @importFrom sf st_cast st_as_sf st_as_sfc st_intersection st_union st_sample st_join
#' st_buffer as_Spatial st_centroid st_set_precision st_make_valid st_set_agr
#' st_geometry st_area st_is_empty st_set_crs st_crs sf_use_s2 st_geometry<-
#' st_drop_geometry st_geometry_type
#'
#' @importFrom dplyr mutate row_number select as_tibble left_join if_else filter
#' arrange desc distinct
#'
#' @importFrom raster raster extend extent focal res crs mask crop rasterize
#' rasterToPoints rasterToPolygons rasterFromXYZ aggregate values ncell
#' values<- adjacent resample ncol
#'
#' @importFrom dplyr select filter mutate
#' @importFrom tidyr unnest
#' @importFrom tibble add_column
#' @importFrom sp proj4string<- coordinates<-
#' @importFrom sf st_as_sf st_union st_distance st_crs
#' @importFrom lwgeom st_snap_to_grid
#' @importFrom smoothr smooth
#' @importFrom gdistance transition geoCorrection
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @importFrom stats na.exclude
#'
#' @export
#'
#' @examples
#' set.seed(1)
#' data(Paracou6_2016)
#' data(DTMParacou)
#' data(PlotMask)
#' data(SpeciesCriteria)
#' data(HarvestableAreaOutputsCable)
#' data(MainTrails)
#'
#' scenario <- "manual"
#' winching <- "2"
#' objective <- 10
#' fuel <- "2"
#' diversification <- TRUE
#'
#' inventory <- addtreedim(cleaninventory(Paracou6_2016, PlotMask),
#' volumeparameters = ForestZoneVolumeParametersTable)
#'
#' treeselectionoutputs <- suppressWarnings(treeselection(inventory,
#' topography = DTMParacou,
#' speciescriteria = SpeciesCriteria,
#' scenario = "manual", objective = objective,
#' fuel = fuel,
#' diversification = diversification,
#' winching = winching,
#' specieslax = FALSE, objectivelax = TRUE,
#' harvestablearea = HarvestableAreaOutputsCable$HarvestableArea,
#' plotslope = HarvestableAreaOutputsCable$PlotSlope,
#' maintrails = MainTrails,
#' harvestablepolygons = HarvestableAreaOutputsCable$HarvestablePolygons,
#' advancedloggingparameters = loggingparameters()))
#'
#' \dontrun{
#' secondtrails <- secondtrailsopening(
#' topography = DTMParacou,
#' plotmask = PlotMask,
#' maintrails = MainTrails,
#' plotslope = HarvestableAreaOutputsCable$PlotSlope,
#' harvestablepolygons = HarvestableAreaOutputsCable$HarvestablePolygons,
#' machinepolygons = HarvestableAreaOutputsCable$MachinePolygons,
#' treeselectionoutputs = treeselectionoutputs,
#' scenario = "manual",
#' winching = winching,
#' fuel = "0",
#' advancedloggingparameters = loggingparameters())
#' }
#'
#' data(SecondaryTrails)
#' library(ggplot2)
#' library(sf)
#'
#' NewInventory <- SecondaryTrails$inventory
#'
#'
#' Harvestable <- sf::st_as_sf(
#' dplyr::filter(NewInventory, LoggingStatus == "harvestable"),
#' coords = c("Xutm", "Yutm")) %>%
#' st_set_crs(st_crs(SecondaryTrails$MainTrailsAccess))
#'
#' HarvestableUp <- sf::st_as_sf(
#' dplyr::filter(NewInventory, LoggingStatus == "harvestableUp"),
#' coords = c("Xutm", "Yutm")) %>%
#' st_set_crs(st_crs(SecondaryTrails$MainTrailsAccess))
#'
#' Selected <- sf::st_as_sf(
#' dplyr::filter(NewInventory, Selected == "1"), coords = c("Xutm", "Yutm")) %>%
#' st_set_crs(st_crs(SecondaryTrails$MainTrailsAccess))
#'
#' Reserve <- sf::st_as_sf(
#' dplyr::filter(NewInventory, LoggingStatus == "reserve"),
#' coords = c("Xutm", "Yutm")) %>%
#' st_set_crs(st_crs(SecondaryTrails$MainTrailsAccess))
#'
#' Future <- sf::st_as_sf(
#' dplyr::filter(NewInventory, LoggingStatus == "future"),
#' coords = c("Xutm", "Yutm")) %>%
#' st_set_crs(st_crs(SecondaryTrails$MainTrailsAccess))
#'
#' ProbedHollow <- sf::st_as_sf(
#' dplyr::filter(NewInventory, ProbedHollow == "1"),
#' coords = c("Xutm", "Yutm")) %>%
#' st_set_crs(st_crs(SecondaryTrails$MainTrailsAccess))
#'
#' VisibleDefect <- sf::st_as_sf(
#' dplyr::filter(NewInventory, VisibleDefect == "1"),
#' coords = c("Xutm", "Yutm")) %>%
#' st_set_crs(st_crs(SecondaryTrails$MainTrailsAccess))
#'
#'
#' ggplot() +
#' # Harvestable zones
#' geom_sf(data = HarvestableAreaOutputsCable$HarvestablePolygons,
#' fill = "olivedrab", alpha = 0.1) +
#' geom_sf(data = HarvestableAreaOutputsCable$MachinePolygons,
#' fill = "olivedrab", alpha = 0.5) +
#' labs(alpha = "Harvestable") +
#' labs(title = "Paracou P6 - secondary trails") +
#'
#' geom_sf(data = VisibleDefect,
#' aes(colour = "Visible defect"), show.legend = "point") +
#' # geom_sf(data = Future,
#' # aes(colour = "Future"), show.legend = "point", size = 4) +
#' geom_sf(data = Reserve,
#' aes(colour = "Reserve"), show.legend = "point", size = 4) +
#' geom_sf(data = Harvestable,
#' aes(colour = "Harvestable"), show.legend = "point", size = 4) +
#' geom_sf(data = HarvestableUp,
#' aes(colour = "HarvestableUp"), show.legend = "point", size = 4) +
#' geom_sf(data = Selected,
#' aes(colour = "Selected"), show.legend = "point") +
#' geom_sf(data = ProbedHollow,
#' aes(colour = "Probed hollow"), show.legend = "point") +
#'
#' # 2ndary trails
#' geom_sf(data = st_as_sf(SecondaryTrails$SmoothedTrails), col = "darkgreen") +
#' geom_sf(data = st_as_sf(SecondaryTrails$MainTrailsAccess), col = "black") +
#'
#' scale_colour_manual(values = c("Visible defect" = "pink",
#' "Harvestable" = "skyblue", "HarvestableUp" = "blue", "Selected" = "red",
#' "Future" = "orange", "Reserve" = "purple", "Probed hollow" = "forestgreen",
#' "Second trails" = "darkgreen", "Harvestable area" = "olivedrab"))
#'
#' SecondaryTrails$TrailsIdentity
#'
secondtrailsopening <- function(
topography,
plotmask,
maintrails,
plotslope,
harvestablepolygons,
machinepolygons,
maintrailsaccess = NULL,
treeselectionoutputs,
scenario,
winching = NULL,
fuel = NULL,
verbose = FALSE,
advancedloggingparameters = loggingparameters()
){
# Arguments check
if(!inherits(treeselectionoutputs, "list"))
stop("The 'treeselectionoutputs' arguments of the 'secondtrailsopening'
function must be list following treeselection output format")
if(!inherits(plotmask, "SpatialPolygonsDataFrame"))
stop("The 'plotmask' argument of the 'secondtrailsopening' function must be
SpatialPolygonsDataFrame")
# if(!any(unlist(lapply(list(maintrails), inherits, "sf" ))))
# stop("The 'maintrails' argument of the 'secondtrailsopening' function must be sf polygon")
if(!inherits(topography, "RasterLayer"))
stop("The 'topography' argument of the 'secondtrailsopening' function must
be RasterLayer")
if(sf::st_is_empty(treeselectionoutputs$SelectedTreesPoints)[1]) {
stop("The 'treeselectionoutputs' argument does not contain any selected
tree.")
}
# Options
options("rgdal_show_exportToProj4_warnings"="none") # to avoid gdal warnings
# gc() # remove intermediary files (to avoid gdal warnings)
# Global Variables
slope <- x <- y <- Harvestable <- idTree <- ID <- type <- ptAcc <- . <- NULL
EstCost <- n.overlaps <- TypeAcc <- IDpts <- Logged <- AccessPolygons <- NULL
Selected <- DeathCause <- ID_Acc <- isEmpty <- gprlAcc <- cblAcc <- NULL
ID.y <- IdMachineZone <- IdMachineZone.y <- IdMachineZone.x <- LineID <- LoggedTrees <- TypeExpl <- NULL
#### Redefinition of the parameters according to the chosen scenario ####
scenariosparameters <- scenariosparameters(scenario = scenario, winching = winching, fuel = fuel)
WinchingInit <- scenariosparameters$winching
winching <- WinchingInit
fuel <- scenariosparameters$fuel
##################################
suppressMessages(sf_use_s2(FALSE)) # to deal with actual unresolved s2 issues in sf ("Spherical geometry (s2) switched off")
CostMatrix <- advancedloggingparameters$CostMatrix
factagg <- floor(advancedloggingparameters$SlopeDistance/res(topography)[1])
advancedloggingparameters$CableLength <- advancedloggingparameters$CableLength + factagg
# Transformation of the DTM so that the maintrails are outside the plot
DTMExtExtended <- raster::extend(topography, c(2*factagg,2*factagg)) # extend the extent
fill.boundaries <- function(x) {
center = 0.5 + (factagg*factagg/2)
if( is.na(x)[center] ) {
return( round(mean(x, na.rm=T),5) )
} else {
return( x[center] )
}
}
DTMExtended <- raster::focal(DTMExtExtended,
matrix(1,factagg,factagg),
fun=fill.boundaries,
na.rm=F, pad=T)
for (iterExt in 1:5) {
DTMExtended <- raster::focal(DTMExtended,
matrix(1,factagg,factagg),
fun=fill.boundaries,
na.rm=F, pad=T)
}
# Generate accessible area from HarvestablePolygones and winching choice
AccessPolygons <- machinepolygons
maintrailsRed <- maintrails %>% st_as_sfc() %>%
st_cast("POLYGON") %>% st_buffer(-factagg) %>%
st_cast("LINESTRING") %>% st_buffer(factagg+0.5)
#-------------------------------------------
if (!is.null(maintrailsaccess)) {
AccessMainTrails <- AccessPolygons %>% st_union() %>%
st_cast("POLYGON", warn = FALSE) %>%
st_as_sf() %>% st_join(maintrailsaccess) %>%
select(ID)
PartMainTrails <- st_intersection(st_geometry(maintrailsRed),
st_geometry(AccessMainTrails)) %>%
st_union(by_feature = T) %>% st_as_sf() %>%
filter(st_geometry_type(.) %in% c("POLYGON","MULTIPOLYGON") ) %>%
st_cast("MULTIPOLYGON", warn = FALSE) %>% # "In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only"
st_as_sf() %>%
st_set_agr(value = "constant") %>%
st_join(AccessMainTrails)%>%
st_make_valid()
PartMainTrails <- PartMainTrails %>%
arrange(desc(st_area(PartMainTrails))) %>%
filter(duplicated(PartMainTrails$ID) == FALSE)
AccessPointAll <- maintrailsaccess
}else{
accesspts <- genaccesspts(topography = topography,
machinepolygons = machinepolygons,
maintrails = maintrails,
advancedloggingparameters = advancedloggingparameters)
PartMainTrails <- accesspts$PartMainTrails
# Generate point access in the intersections between accessible area and maintrails (ID = accessible area index)
AccessPointAll <- accesspts$AccessPointAll
AccessMainTrails <- AccessPolygons %>% st_union() %>%
st_cast("POLYGON", warn = FALSE) %>%
st_as_sf() %>% st_join(accesspts$AccessPointAll) %>%
select(ID)
}
#-------------------------------------------
# Generate Cost raster --> cf CostMatrix
CostRaster <- raster(extent(DTMExtended),resolution = res(DTMExtended), crs = crs(DTMExtended))
values(CostRaster) <- CostMatrix[[2]][[1]]$CostValue
#Generate weight according to slope conditions
# RastersToPoints
plotslopePoint <- as_tibble(rasterToPoints(plotslope))
CostRasterPoint <- as_tibble(rasterToPoints(CostRaster))
# left join par x et y
plotTib <-
left_join(plotslopePoint, CostRasterPoint, by = c('x', 'y'))
CostSlopeRaster <- plotTib %>%
mutate(Harvestable = if_else(
slope <= atan(CostMatrix[[1]][[1]]$Slope/100),
true = CostMatrix[[1]][[1]]$Cost,
false = if_else(
slope <= atan(CostMatrix[[1]][[2]]$Slope/100),
true = CostMatrix[[1]][[2]]$Cost,
false = if_else(
slope <= atan(CostMatrix[[1]][[3]]$Slope/100),
true = CostMatrix[[1]][[3]]$Cost,
false = if_else(
slope <= atan(CostMatrix[[1]][[4]]$Slope/100),
true = CostMatrix[[1]][[4]]$Cost,
false = if_else(
slope <= atan(CostMatrix[[1]][[5]]$Slope/100),
true = CostMatrix[[1]][[5]]$Cost,
false = CostMatrix[[1]][[6]]$Cost)
)
)
)
)) %>% dplyr::select(x,y,Harvestable)
CostRaster <- rasterFromXYZ(CostSlopeRaster, crs = crs(DTMExtended)) # Update weights from plotTib tibble
#Aggregation each raster to selected resolution
CostRasterMean <- aggregate(CostRaster, fact=factagg, fun=mean)
DTMmean <- aggregate(DTMExtended, fact=factagg, fun=mean)
plotmaskagg <- rasterToPolygons(DTMmean, dissolve=T) %>% st_as_sf() %>% st_union() %>% as_Spatial()
#Generate protection buffer for Big Trees (dist : ScndTrailWidth/2 + 2 m)
BigTreesRaster <- raster(extent(DTMmean),resolution = res(DTMmean), crs = crs(DTMmean))
values(BigTreesRaster) <- 0
BigTreesRaster <- rasterize(x = as_Spatial(treeselectionoutputs$BigTreesPoints %>%
st_buffer(dist = advancedloggingparameters$ScndTrailWidth/2 + 2) ),
y = BigTreesRaster,
field = CostMatrix[[2]][[3]]$CostValue,
update = TRUE)
BigTreesRaster <- crop(BigTreesRaster, CostRasterMean)
BigTreesRaster <- mask(BigTreesRaster, plotmaskagg)
#Generate protection buffer for Reserve Trees (dist : ScndTrailWidth/2 + 2 m)
ReserveRaster <- raster(extent(DTMmean),resolution = res(DTMmean), crs = crs(DTMmean))
values(ReserveRaster) <- 0
ReserveRaster <- rasterize(x = as_Spatial(treeselectionoutputs$HarvestableTreesPoints %>%
st_buffer(dist = advancedloggingparameters$ScndTrailWidth/2 + 2) ),
y = ReserveRaster,
field = CostMatrix[[2]][[4]]$CostValue,
update = TRUE)
ReserveRaster <- crop(ReserveRaster, CostRasterMean)
ReserveRaster <- mask(ReserveRaster, plotmaskagg)
#Generate protection buffer for Futures Trees (dist : ScndTrailWidth/2 + 2 m)
FutureRaster <- raster(extent(DTMmean),resolution = res(DTMmean), crs = crs(DTMmean))
values(FutureRaster ) <- 0
FutureRaster <- rasterize(x = as_Spatial(treeselectionoutputs$HarvestableTreesPoints %>%
st_buffer(dist = advancedloggingparameters$ScndTrailWidth/2 +2) ),
y = FutureRaster ,
field = CostMatrix[[2]][[5]]$CostValue,
update = TRUE)
FutureRaster <- crop(FutureRaster, CostRasterMean)
FutureRaster <- mask(FutureRaster, plotmaskagg)
#Update Cost Raster with protection buffers
#Generate protection buffer for selected Trees (dist : ScndTrailWidth/2 + 2 m)
SelectedRaster <- raster(extent(DTMmean),resolution = res(DTMmean), crs = crs(DTMmean))
values(SelectedRaster ) <- 0
SelectedRaster <- rasterize(x = as_Spatial(treeselectionoutputs$SelectedTreesPoints %>%
st_buffer(dist = 2) ),
y = SelectedRaster,
field = CostMatrix[[2]][[3]]$CostValue/2,
update = TRUE)
SelectedRaster <- crop(SelectedRaster, CostRasterMean)
SelectedRaster <- mask(SelectedRaster, plotmaskagg)
#Update Cost Raster with protection buffers
CostRasterMean <- CostRasterMean + BigTreesRaster + ReserveRaster + FutureRaster + SelectedRaster
#Generate maintrail intersection cost raster
CostRasterMean <- rasterize(x = as_Spatial(AccessPointAll %>% st_buffer(dist = advancedloggingparameters$ScndTrailWidth+2)),
y = CostRasterMean ,
field = CostMatrix[[2]][[6]]$CostValue,
update = TRUE)
#Generate Slope accessibility for grapple machine
if (winching == "2") {
plotslopePointGrpl <- as_tibble(rasterToPoints(plotslope))
CostRasterPointGrpl <- as_tibble(rasterToPoints(CostRaster))
# left join par x et y
plotTibGrpl <-
left_join(plotslopePointGrpl, CostRasterPointGrpl, by = c('x', 'y'))
CostSlopeRasterGrpl <- plotTibGrpl %>%
mutate(Harvestable = if_else(
slope <= atan(advancedloggingparameters$GrappleMaxslope/100),
true = 0,
false = Inf
)) %>% dplyr::select(x,y,Harvestable)
CostRasterGrpl <- rasterFromXYZ(CostSlopeRasterGrpl, crs = crs(CostRasterMean))
PolygonGrpl <- rasterToPolygons(x = CostRasterGrpl,
n = 8,
dissolve = TRUE)
CostRasterMeanGrpl <- rasterize(x = PolygonGrpl[PolygonGrpl$Harvestable == Inf,],
y = CostRasterMean ,
field = Inf,
update = TRUE)
CostRasterMeanGrpl <- rasterize(x = as_Spatial(AccessPointAll %>% st_buffer(dist = advancedloggingparameters$ScndTrailWidth+2)),
y = CostRasterMeanGrpl ,
field = CostMatrix[[2]][[6]]$CostValue,
update = TRUE)
#Generate accessible weights raster
AccessRaster <- raster(extent(CostRaster),resolution = res(CostRaster), crs = crs(CostRaster))
values(AccessRaster) <- CostMatrix[[2]][[2]]$CostValue
AccessRaster <- rasterize(x = as_Spatial(machinepolygons %>% st_buffer(dist = factagg)),
y = AccessRaster ,
field = 1,
update = TRUE)
AccessRaster <- aggregate(AccessRaster,fact=factagg, fun=min)
AccessRaster <- crop(AccessRaster, CostRasterMeanGrpl)
AccessRaster <- mask(AccessRaster, plotmaskagg)
#Update Cost Raster with accessible weights raster
CostRasterMeanGrpl <- CostRasterMeanGrpl * AccessRaster
}
#Generate accessible weights raster
AccessRaster <- raster(extent(CostRaster),resolution = res(CostRaster), crs = crs(CostRaster))
values(AccessRaster) <- CostMatrix[[2]][[2]]$CostValue
AccessRaster <- rasterize(x = as_Spatial(machinepolygons %>% st_buffer(dist = factagg)),
y = AccessRaster ,
field = 1,
update = TRUE)
AccessRaster <- aggregate(AccessRaster,fact=factagg, fun=min)
AccessRaster <- crop(AccessRaster, CostRasterMean)
AccessRaster <- mask(AccessRaster,plotmaskagg)
#Update Cost Raster with accessible weights raster
CostRasterMean <- CostRasterMean * AccessRaster
pathLines <- list() #Initialize storage pathlines
Lines <- list() #Initialize storage logged trees
k <- 1 #Initialize pathlines index
l <- 1 #Initialize lines index
#Generate appropriate selected trees points format
ptsAllinit <- treeselectionoutputs$SelectedTreesPoints %>%
dplyr::select(idTree) %>%
st_cast("POINT", warn = FALSE) %>%
mutate(ID = NA) %>%
mutate(type = "Tree") %>%
st_set_crs(st_crs(AccessPointAll)) %>%
st_join(AccessMainTrails %>% st_buffer(dist = max((advancedloggingparameters$ScndTrailWidth/2),factagg)) %>%
mutate(ID_Acc = ID) %>% dplyr::select(-ID)) %>%
dplyr::select(ID,ID_Acc,type,idTree)
ptsAllNA <- ptsAllinit %>%
filter(is.na(ID_Acc)) %>%
mutate(ID_Acc = AccessMainTrails$ID[c(max.col(-st_distance(ptsAllinit %>%
filter(is.na(ID_Acc)),AccessMainTrails)))])
ptsAllNotNA <- ptsAllinit %>%
filter(!is.na(ID_Acc))
ptsAll <- rbind(ptsAllNotNA,ptsAllNA)
# Reassign Selected Tree values (= BigTrees) to the aggregated Cost raster
CostRasterMean <- rasterize(x = as_Spatial(ptsAll %>% st_buffer(dist = factagg)),
y = CostRasterMean ,
field = CostMatrix[[2]][[3]]$CostValue,
update = TRUE)
if (winching == "2") {
# Reassign Selected Tree values (= BigTrees) to the aggregated Cost raster (Grpl)
CostRasterMeanGrpl <- rasterize(x = as_Spatial(ptsAll %>%
st_buffer(dist = factagg)),
y = CostRasterMeanGrpl ,
field = CostMatrix[[2]][[3]]$CostValue,
update = TRUE)
CostRasterMeanGrpl <- crop(CostRasterMeanGrpl, CostRasterMean)
CostRasterMeanGrpl <- mask(CostRasterMeanGrpl, plotmaskagg)
}
#Generate maintrail intersection cost raster
CostRasterMean <- rasterize(x = as_Spatial(AccessPointAll %>% st_buffer(dist = 3*factagg)),
y = CostRasterMean ,
field = CostMatrix[[2]][[6]]$CostValue,
update = TRUE)
CostRasterMean <- crop(CostRasterMean,DTMmean)
CostRasterMean <- mask(CostRasterMean,plotmaskagg)
DTMmean <- crop(DTMmean,CostRasterMean)
DTMmean <- mask(DTMmean,plotmaskagg)
#Compute conductance raster
CondSurf <- 1/CostRasterMean
########### Compute slope criteria transition graph ###############
if (winching == "2") {
#Compute adjacent transition layer according to slope conditions (winching = "2")
SlopeCondGrpl <- sloperdcond(topography = DTMmean,
advancedloggingparameters = advancedloggingparameters,
grapple = TRUE)
}
#Compute adjacent transition layer according to slope conditions (winching = "1")
SlopeCond <- sloperdcond(topography = DTMmean,advancedloggingparameters = advancedloggingparameters)
########### Compute LCP algorithm ###############
for(ID_Access in unique(ptsAll$ID_Acc)) { # ici
winching <- WinchingInit
pts <- ptsAll %>%
filter(ID_Acc == ID_Access) %>% dplyr::mutate(ID = ID_Access) %>%
dplyr::select(-ID_Acc) %>%
dplyr::distinct()
AccessPoint <- AccessPointAll %>%
filter(ID == ID_Access)
pathLinesWIP <- c()
ki <- 1
if (winching == "0") {
# winching 0 #########
#Update Cost Raster with accessible weights raster
CostRasterMean <- CostRasterMean * AccessRaster
CondSurf <- 1/CostRasterMean
pts <- st_set_crs(pts, st_crs(AccessPoint))# set crs
pts <- rbind(AccessPoint %>%
mutate(type = "Access") %>%
mutate(IPpts = paste0("A.",row_number())),
pts %>%
mutate(IPpts = paste0("T.",idTree)))
PointAcc <- pts %>%
filter(type == "Access") %>%
mutate(EstCost = NA)
#Compute Cost between points and Access points
CostDistEst <- adjtopolcp(costSurface = CondSurf,
topography = DTMmean ,
pts = as_Spatial(pts),
slopeRdCond = SlopeCond,
paths = FALSE) [,1:dim(AccessPoint)[1]]
CostDistEst <- CostDistEst[(length(PointAcc$ID)+1):length(CostDistEst)[1]]
CostDistEst <- matrix(CostDistEst,ncol = length(PointAcc$ID) )
#Attribute a least cost point access to each points
PointTree <- pts %>% filter(type == "Tree") %>%
mutate(ptAcc = max.col(-CostDistEst))
PointTree$ID <- as.factor(PointTree$ptAcc)
levels(PointTree$ID) <- as.character(AccessPoint$ID)
pts <- rbind(AccessPoint %>%
mutate(type = "Access") %>%
mutate(ptAcc = row_number()) %>%
mutate(IPpts = paste0("A.",row_number())),
PointTree %>%
mutate(IPpts = paste0("T.",idTree)))
PointAcc <- pts %>%
filter(type == "Access") %>%
mutate(EstCost = NA)
#Attribute a least cost point access to each points
PointTree <- pts %>% filter(type == "Tree") %>%
mutate(ptAcc = max.col(-CostDistEst))
if (length(PointAcc$ID) == 1) {
PointTree <- PointTree %>% mutate(EstCost = CostDistEst)
}else{PointTree <- PointTree %>% mutate(EstCost = apply(CostDistEst,1, min))}
selectedPoints <- rbind(PointAcc,PointTree)
accessPtId <- 1
TmpSelectedPts <- selectedPoints %>%
filter(ptAcc == accessPtId)
TmpAccPts <- TmpSelectedPts %>%
filter(type == "Access")
#Select Costliest selected tree
TmpTreePts <- TmpSelectedPts %>%
filter(type == "Tree") %>%
arrange(desc(EstCost))
if (dim(TmpTreePts)[1] != 0) {
for (TreeId in 1:dim(TmpTreePts)[1]){
#Compute Least cost path polygons to the WIP selected tree
TmpPtsWIP <- rbind(TmpAccPts,TmpTreePts[TreeId,])
TmpPathWIP <- adjtopolcp(costSurface = CondSurf,topography = DTMmean , pts = TmpPtsWIP %>% as_Spatial(),
slopeRdCond = SlopeCond,paths = TRUE)
if (TmpPathWIP[[1]][1,2] == 0) {
#Update Cost raster with LCP second trail
CostRasterMean <- rasterize(x = TmpTreePts[TreeId,] ,
y = CostRasterMean ,
field = CostMatrix[[2]][[7]]$CostValue,update = TRUE)
}else{
#Update Cost raster with LCP second trail
CostRasterMean <- rasterize(x = TmpPathWIP[[2]] ,
y = CostRasterMean ,
field = CostMatrix[[2]][[7]]$CostValue,update = TRUE)
}
CondSurf <- 1/CostRasterMean
#Store pathline
pathLines[[k]] <- TmpPathWIP[[2]]
pathLines[[k]]@lines[[1]]@ID <- paste("Path", TmpPtsWIP$idTree[2], sep = ".")
Lines[[k]] <- list("LineID" = k,"LoggedTrees" = TmpPtsWIP$idTree[2],"TypeExpl" = "FoT")
k <- k +1
}
}
}else{
# winching 1/2 #########
# Select intersection points from buffer polygons
PointAcc <- AccessPoint %>% #def Access Point
mutate(type = "Access") %>%
mutate(IDpts = paste0("A.",row_number())) %>%
mutate(n.overlaps = NA, origins = NA) %>%
dplyr::select(-ID)
TreePts <- pts %>% filter(type == "Tree")
if (winching == "2") {
ptsGrpl <- TreePts %>% #def Grpl point
st_buffer(dist = advancedloggingparameters$GrappleLength) %>%
st_snap_to_grid(size = .2) %>% # avoid GEOS error (st-intersection issue)
#st_set_precision(1) %>%
st_temporaryintersection(topography = topography,
plotmask = plotmask,
advancedloggingparameters = advancedloggingparameters) %>%
st_make_valid() %>%
mutate(type = "Inter") %>%
dplyr::select(-ID)
for (OriginI in 1:length(ptsGrpl$origins)) {
for (OriginJ in 1:length(ptsGrpl$origins[[OriginI]])) {
IndexPts <- ptsGrpl$origins[[OriginI]][OriginJ]
ptsGrpl$origins[[OriginI]][OriginJ] <- TreePts$idTree[IndexPts]
}
}
ptsGrpl$isEmpty <- FALSE
for (h in 1:dim(ptsGrpl)[1]) {
if (ptsGrpl[h,] %>% st_intersects(AccessMainTrails %>%
filter(ID == ID_Access),sparse = FALSE) == FALSE) {
ptsGrpl$isEmpty[h] <- TRUE
}else{
suppressWarnings(st_geometry(ptsGrpl[h,]) <- st_geometry(ptsGrpl[h,] %>% st_intersection(AccessMainTrails %>%
filter(ID == ID_Access))))
}
}
#Filter polygons which intersect accessible area to second trails
ptsGrpl<- st_set_crs(ptsGrpl,st_crs(AccessPolygons)) # set crs from AccessPolygons
ptsGrpl <- ptsGrpl %>%
filter(isEmpty == FALSE) %>%
dplyr::select(-isEmpty) %>%
mutate(IDpts = paste0("I.", row_number()))
ptsWIP <- ptsGrpl %>% #def Grpl point as WIP point
st_set_agr("constant") %>%
st_centroid()
ptsCbl <- TreePts %>% #def cbl polygons
st_buffer(dist = advancedloggingparameters$CableLength) %>%
st_snap_to_grid(size = 1) %>%# avoid GEOS error (st-intersection issue)
st_set_precision(1) %>%
st_temporaryintersection(topography = topography,
plotmask = plotmask,
advancedloggingparameters = advancedloggingparameters) %>%
st_make_valid() %>%
mutate(type = "Inter")%>%
dplyr::select(-ID)
for (OriginI in 1:length(ptsCbl$origins)) {
for (OriginJ in 1:length(ptsCbl$origins[[OriginI]])) {
IndexPts <- ptsCbl$origins[[OriginI]][OriginJ]
ptsCbl$origins[[OriginI]][OriginJ] <- TreePts$idTree[IndexPts]
}
}
ptsCbl$isEmpty <- FALSE
for (h in 1:dim(ptsCbl)[1]) {
if (ptsCbl[h,] %>% st_intersects(AccessMainTrails %>%
filter(ID == ID_Access),sparse = FALSE) == FALSE) {
ptsCbl$isEmpty[h] <- TRUE
}else{
suppressWarnings(st_geometry(ptsCbl[h,]) <- st_geometry(ptsCbl[h,] %>%
st_intersection(AccessMainTrails %>%
filter(ID == ID_Access))))
}
}
ptsCbl <- st_set_crs(ptsCbl,st_crs(AccessPolygons)) #set crs from AccessPolygons
ptsCbl <- ptsCbl %>% filter(isEmpty == FALSE) %>% dplyr::select(-isEmpty) %>%
mutate(IDpts = paste0("I.",row_number())) #Filter polygons which intersect accessible area to second trails
ptsWIPCbl <- ptsCbl %>%#Convert polygons to centroid
st_set_agr("constant") %>%
st_centroid()
ptsWIP <- ptsWIP %>%
arrange(desc(n.overlaps))
RemainTree <- max(length(unique(unlist(ptsWIP$origins))),length(unique(unlist(ptsWIPCbl$origins))))
#Select adjacent grpl graph
SlopeCondRd <- SlopeCondGrpl
}else{
ptsCbl <- TreePts %>% #def cbl point
st_buffer(dist = advancedloggingparameters$CableLength) %>%
st_snap_to_grid(size = 1) %>% # avoid GEOS error (st-intersection issue)
st_set_precision(1) %>%
st_temporaryintersection(topography = topography,
plotmask = plotmask,
advancedloggingparameters = advancedloggingparameters) %>%
st_make_valid() %>%
mutate(IDpts = paste0("I.",row_number())) %>%
mutate(type = "Inter") %>% dplyr::select(-ID)
ptsCbl <- st_set_crs(ptsCbl,st_crs(AccessPolygons)) #set crs from AccessPolygons
for (OriginI in 1:length(ptsCbl$origins)) {
for (OriginJ in 1:length(ptsCbl$origins[[OriginI]])) {
IndexPts <- ptsCbl$origins[[OriginI]][OriginJ]
ptsCbl$origins[[OriginI]][OriginJ] <- TreePts$idTree[IndexPts]
}
}
ptsCbl$isEmpty <- FALSE
for (h in 1:dim(ptsCbl)[1]) {
if (st_geometry(ptsCbl[h,]) %>% st_intersects(AccessMainTrails %>%
filter(ID == ID_Access),sparse = FALSE) == FALSE) {
ptsCbl$isEmpty[h] <- TRUE
}else{
suppressWarnings(st_geometry(ptsCbl[h,]) <- st_geometry(ptsCbl[h,] %>%
st_intersection(AccessMainTrails %>%
filter(ID == ID_Access))))
ptsCbl$isEmpty[h] <- FALSE
}
}
ptsWIP <- ptsCbl %>% filter(isEmpty == FALSE) %>%
dplyr::select(-isEmpty) %>% #def cbl point as WIP point
st_set_agr("constant") %>%
st_centroid()
ptsWIP <- ptsWIP %>% #filter cbl intersection centroid point out plot
filter(st_intersects(st_geometry(ptsWIP),
st_geometry(plotmask %>% st_as_sf()),
sparse = FALSE)) %>%
arrange(desc(n.overlaps))
RemainTree <- length(unique(ptsWIP$idTree))
#Select adjacent cbl graph
SlopeCondRd <- SlopeCond
}
#Define a switch from grapple accessible tree exploitation to cable exploitation
if (winching == "2") {
WinchingSwitch <- FALSE
}else{
WinchingSwitch <- TRUE
}
steps <- 0
RemainTreeInit <- RemainTree
#Loop over possible intersection
while (RemainTree !=0L & steps < 5*RemainTreeInit) {
steps <- steps +1
if (verbose) {
print(RemainTree)
}
#Switch from grpl to cbl exploitation when grapple accessible tree != 0
Grpl2CblFlag <- FALSE
if (length(unique(unlist(ptsWIP$origins))) == 0 & WinchingSwitch == FALSE & winching == "2") {
winching <- "1"
ptsWIP <- ptsWIPCbl
WinchingSwitch <- TRUE
}
ptsWIPmax <- rbind(PointAcc,ptsWIP %>%
filter(st_intersects(.,st_geometry(plotmask %>%
st_as_sf()),sparse = FALSE)) %>%
filter(n.overlaps == max(ptsWIP$n.overlaps))) %>%
mutate(TypeAcc = NA) %>%
mutate(EstCost = NA) %>%
mutate(ptsAcc = NA)
ptsTreeWIP <- rbind(PointAcc,TreePts %>%
mutate(n.overlaps = NA, origins = idTree,IDpts = NA) %>%
dplyr::select(-ID)) %>%
mutate(TypeAcc = NA) %>%
mutate(EstCost = NA) %>%
mutate(ptsAcc = NA)
if (ki == 1 ) {
if (winching == "2") {
ptsDirAcc <- ptsTreeWIP %>%
mutate(gprlAcc = c(FALSE,as.numeric(st_distance(ptsTreeWIP)[2:dim(ptsTreeWIP)[1],1]) < advancedloggingparameters$GrappleLength)) %>%
filter(gprlAcc == TRUE | type == "Access") %>%
dplyr::select(-gprlAcc)
TmpTypeAcc <- "Grpl"
}else{
ptsDirAcc <- ptsTreeWIP %>%
mutate(cblAcc = c(FALSE,as.numeric(st_distance(ptsTreeWIP)[2:dim(ptsTreeWIP)[1],1]) < advancedloggingparameters$CableLength)) %>%
filter(cblAcc == TRUE | type == "Access") %>%
dplyr::select(-cblAcc)
TmpTypeAcc <- "Cbl"
}
}else{
# Set pathsWIP CRS
pathsWIP <- st_set_crs(pathsWIP, st_crs(PointAcc))
if (winching == "2") {
ptsDirAcc <- ptsTreeWIP %>%
mutate(gprlAcc = c(FALSE,as.numeric(st_distance(ptsTreeWIP,st_union(pathsWIP %>%
st_buffer(dist = advancedloggingparameters$MaxMainTrailWidth/2) %>%
st_make_valid() %>%
st_union(),PointAcc %>%
st_buffer(dist = advancedloggingparameters$MaxMainTrailWidth/2) %>%
st_make_valid() %>%
st_union()))[2:dim(ptsTreeWIP)[1],1]) < advancedloggingparameters$GrappleLength)) %>%
filter(gprlAcc == TRUE | type == "Access") %>%
dplyr::select(-gprlAcc)
TmpTypeAcc <- "Grpl"
}else{
ptsDirAcc <- ptsTreeWIP %>% mutate(cblAcc = c(FALSE,as.numeric(st_distance(ptsTreeWIP,st_union(pathsWIP %>% st_buffer(dist = advancedloggingparameters$MaxMainTrailWidth/2) %>% st_make_valid() %>%
st_union(),PointAcc %>% st_buffer(dist = advancedloggingparameters$MaxMainTrailWidth/2) %>% st_make_valid() %>%
st_union()))[2:dim(ptsTreeWIP)[1],1]) < advancedloggingparameters$CableLength)) %>%
filter(cblAcc == TRUE | type == "Access") %>% dplyr::select(-cblAcc)
TmpTypeAcc <- "Cbl"
}
}
if (dim(ptsDirAcc)[1] > 1) {
PointTreeWIP <- ptsDirAcc %>% filter(type == "Tree")
if (verbose){
print(paste0("LoggedTrees : ",PointTreeWIP$origins[[1]]))
}
Lines[[l]] <- list("LineID" = NA,"LoggedTrees" = PointTreeWIP$origins, "TypeExpl" = TmpTypeAcc,"IdMachineZone" = ID_Access)
l <- l+1
}else{
#Compute Cost between points and Access points in cbl exploitation
CostDistEstCbl <- adjtopolcp(costSurface = CondSurf,
topography = DTMmean ,
pts = ptsWIPmax %>%
as_Spatial(),
slopeRdCond = SlopeCond,
paths = FALSE) [,1]
CostDistEstCbl <- CostDistEstCbl[(length(PointAcc$ID)+1):length(CostDistEstCbl)[1]]
CostDistEstCbl <- matrix(CostDistEstCbl,ncol = length(PointAcc$ID) )
#Attribute a least cost point access to each points
PointTree <- ptsWIPmax %>% filter(type != "Access" ) %>%
mutate(ptAccCbl = max.col(- CostDistEstCbl,ties.method = "first")) %>%
mutate(EstCostCbl = CostDistEstCbl)
if (winching == "2") {
CondSurfGrpl <- 1/CostRasterMeanGrpl
#Compute Cost between points and Access points in grpl exploitation
CostDistEstGrpl <- adjtopolcp(costSurface = CondSurfGrpl,
topography = DTMmean ,
pts = ptsWIPmax %>%
as_Spatial(),
slopeRdCond = SlopeCondGrpl,
paths = FALSE) [,1:length(PointAcc$ID)]
CostDistEstGrpl <- CostDistEstGrpl[(length(PointAcc$ID)+1):length(CostDistEstGrpl)[1]]
CostDistEstGrpl <- matrix(CostDistEstGrpl,ncol = length(PointAcc$ID) )
#Attribute a least cost point access to each points
PointTree <- PointTree %>%
mutate(ptAccGpl = max.col(-matrix(CostDistEstGrpl,ncol = length(PointAcc$ID) ),ties.method = "first")) %>%
mutate(EstCostGrpl = CostDistEstGrpl)
#Prioritize grpl exploitation if possible
CostDistEstGrpl[CostDistEstGrpl != Inf] <- 0
CostDistEstGrpl[CostDistEstGrpl == Inf] <- 1
for (j in 1:length(CostDistEstGrpl)) {
PointTree[j,"TypeAcc"] <- as.character(prod(CostDistEstGrpl[j,]))
}
PointTree$TypeAcc[PointTree$TypeAcc == "0"] <- "Grpl"
PointTree$TypeAcc[PointTree$TypeAcc == "1"] <- "Cbl"
PointTree$ptsAcc[PointTree$TypeAcc == "Grpl"] <- PointTree$ptAccGpl[PointTree$TypeAcc == "Grpl"]
PointTree$ptsAcc[PointTree$TypeAcc != "Grpl"] <- PointTree$ptAccCbl[PointTree$TypeAcc != "Grpl"]
PointTree$EstCost[PointTree$TypeAcc == "Grpl"] <- PointTree$EstCostGrpl[PointTree$TypeAcc == "Grpl"]
PointTree$EstCost[PointTree$TypeAcc != "Grpl"] <- PointTree$EstCostCbl[PointTree$TypeAcc != "Grpl"]
PointTree <- PointTree %>% arrange(desc(TypeAcc),EstCost)
#Define WIP point according to possible exploitation type
if (PointTree$TypeAcc[1] == "Grpl") {
#Reconstruct access points + selected tree in grpl exploitation
TmpPtsWIP <- ptsGrpl %>%
filter(IDpts == PointTree$IDpts[1]) %>%
st_union() %>%
st_cast("POINT", warn = FALSE) %>%
st_as_sf() %>%
mutate(type = "Overlay") %>%
mutate(ptsAcc = PointTree$ptsAcc[1]) %>%
mutate(IDpts = PointTree$IDpts[1]) %>%
mutate(origins = PointTree$origins[1])%>%
mutate(n.overlaps = PointTree$n.overlaps[1])
TmpTypeAcc <- "Grpl"
#Select adjacent grpl graph
SlopeCondRd <- SlopeCondGrpl
}else{
#if the selected point is not a grpl accessible point BUT not all grpl accessible point are done
#Shift to cbl exploitation
Grpl2CblFlag <- TRUE
#Reconstruct access points + selected tree in cbl exploitation
TmpPtsWIP <- ptsCbl %>%
filter(IDpts == PointTree$IDpts[1]) %>%
st_union() %>%
st_cast("POINT", warn = FALSE) %>%
st_as_sf() %>%
mutate(type = "Overlay") %>%
mutate(ptsAcc = PointTree$ptsAcc[1]) %>%
mutate(IDpts = PointTree$IDpts[1]) %>%
mutate(origins = PointTree$origins[1])%>%
mutate(n.overlaps = PointTree$n.overlaps[1])
#Select adjacent cbl graph
SlopeCondRd <- SlopeCond
TmpTypeAcc <- "Cbl"
}
}else{
#Reconstruct access points + selected tree in cbl exploitation
TmpPtsWIP <- ptsCbl %>%
filter(IDpts == PointTree$IDpts[1]) %>%
st_union() %>%
st_cast("POINT", warn = FALSE) %>%
st_as_sf() %>%
mutate(type = "Overlay") %>%
mutate(ptsAcc = PointTree$ptsAcc[1]) %>%
mutate(IDpts = PointTree$IDpts[1]) %>%
mutate(origins = PointTree$origins[1]) %>%
mutate(n.overlaps = PointTree$n.overlaps[1])
#Select adjacent cbl graph
SlopeCondRd <- SlopeCond
TmpTypeAcc <- "Cbl"
}
#filter WIP points in accessible scd trail area
# TmpPtsWIP <- st_set_crs(TmpPtsWIP,st_crs(ptsCbl)) # set crs from ptsCbl
TmpPtsWIP <- TmpPtsWIP %>%
filter(st_intersects(
st_geometry(TmpPtsWIP),st_geometry(AccessPolygons %>%
st_union()),sparse = FALSE))
#Reconstruct access points + selected points
TmpPtsWIP <- rbind( PointAcc %>%
st_union() %>%
st_cast("POINT", warn = FALSE)%>%
st_as_sf() %>%
mutate(type = "Access") %>%
mutate(ptsAcc = NA ) %>%
mutate(IDpts = NA) %>%
mutate(origins = NA)%>%
mutate(n.overlaps = NA),TmpPtsWIP)
#Compute Cost between all points and Access points
CostDistEstWIP <- adjtopolcp(costSurface = CondSurf,topography = DTMmean , pts = TmpPtsWIP %>% as_Spatial(),
slopeRdCond = SlopeCondRd,paths = FALSE)[,1]
CostDistEstWIP <- CostDistEstWIP[(length(PointAcc$ID)+1):length(CostDistEstWIP)[1]]
CostDistEstWIP <- matrix(CostDistEstWIP, ncol = length(PointAcc$ID))
PointTreeWIP <- TmpPtsWIP %>%
filter(type == "Overlay") %>%
mutate(EstCost = CostDistEstWIP) %>%
arrange(EstCost)
TmpPtsWIP <- rbind(TmpPtsWIP %>% filter(type == "Access") %>%
mutate(EstCost = NA),PointTreeWIP[1,])
TmpPathWIP <- adjtopolcp(costSurface = CondSurf,topography = DTMmean , pts = TmpPtsWIP %>% as_Spatial(),
slopeRdCond = SlopeCondRd,paths = FALSE)
TmpPathWIPCost <- TmpPathWIP[1:length(PointAcc$ID),length(PointAcc$ID)+1]
LCPathWIP <- max.col(t(-TmpPathWIPCost))
TmpPtsWIP <- rbind(TmpPtsWIP[LCPathWIP,],PointTreeWIP[1,])
TmpPathWIP <- adjtopolcp(costSurface = CondSurf,topography = DTMmean , pts = TmpPtsWIP %>% as_Spatial(),
slopeRdCond = SlopeCondRd,paths = TRUE)
if (TmpPathWIP[[1]][2,1] != 0) {
if (TmpTypeAcc == "Grpl") {
CostRasterMean <- rasterize(x = TmpPathWIP[[2]] ,
y = CostRasterMean ,
field = CostMatrix[[2]][[7]]$CostValue,update = TRUE)
CostRasterMeanGrpl <- rasterize(x = TmpPathWIP[[2]] ,
y = CostRasterMeanGrpl ,
field = CostMatrix[[2]][[7]]$CostValue,update = TRUE)
}else{
CostRasterMean <- rasterize(x = TmpPathWIP[[2]] ,
y = CostRasterMean ,
field = CostMatrix[[2]][[7]]$CostValue,update = TRUE)
}
}
pathLines[[k]] <- TmpPathWIP[[2]]
pathLines[[k]]@lines[[1]]@ID <- paste("Path",
"A",
LCPathWIP,
"NTree",
length(PointTreeWIP$origins[[1]]),
"T",
paste(as.character(unlist(PointTreeWIP$origins[[1]])),
collapse='-'),
sep = ".")
Lines[[l]] <- list("LineID" = k,"LoggedTrees" = PointTreeWIP$origins[[1]], "TypeExpl" = TmpTypeAcc,"IdMachineZone" = ID_Access)
k <- k +1
l <- l+1
pathLinesWIP[[ki]] <- TmpPathWIP[[2]]
pathLinesWIP[[ki]]@lines[[1]]@ID <- paste("Path",
"A",
LCPathWIP,
"NTree",
length(PointTreeWIP$origins[[1]]),
"T",
paste(as.character(unlist(PointTreeWIP$origins[[1]])),
collapse='-'),
sep = ".")
ki <- ki +1
pathsWIP <- do.call(rbind, pathLinesWIP)
pathsWIP <- pathsWIP %>% st_as_sf() %>% st_make_valid() %>%
filter(!st_is_empty(pathsWIP %>%
st_as_sf())) %>%
filter(st_geometry_type(.) %in% c("LINESTRING"))
if(dim(pathsWIP)[1]>0){
pathsWIP <- pathsWIP$geometry %>% as_Spatial()
if ( length(pathsWIP)==1) {
pathsWIP <- pathsWIP %>% st_as_sf() %>%
st_buffer(dist = advancedloggingparameters$ScndTrailWidth/2) %>% st_make_valid() %>%
st_union()
}else{
pathsWIP <- smoothtrails(paths = pathsWIP ,
plotmask = plotmask,
verbose = FALSE,
advancedloggingparameters = advancedloggingparameters)$SmoothedTrails %>%
st_union()
}
}else{
pathsWIP <- pathsWIP %>% st_union()
}
}
TreePts$Logged<- FALSE
for (j in 1:dim(TreePts)[1]) {
TreePts$Logged[j] <- any(PointTreeWIP$origins[[1]] %in% TreePts$idTree[j])
}
TreePts <- TreePts %>%
filter(Logged == FALSE) %>%
dplyr::select(-Logged)
if (dim(TreePts)[1]> 0) {
if (winching == "2") {
ptsGrpl <- TreePts %>% #def Grpl point
st_buffer(dist = advancedloggingparameters$GrappleLength) %>%
st_snap_to_grid(size = .2) %>% # avoid GEOS error (st-intersection issue)
#st_set_precision(1) %>%
st_temporaryintersection(topography = topography,
plotmask = plotmask,
advancedloggingparameters = advancedloggingparameters) %>%
st_make_valid() %>%
mutate(type = "Inter") %>%
dplyr::select(-ID)
for (OriginI in 1:length(ptsGrpl$origins)) {
for (OriginJ in 1:length(ptsGrpl$origins[[OriginI]])) {
IndexPts <- ptsGrpl$origins[[OriginI]][OriginJ]
ptsGrpl$origins[[OriginI]][OriginJ] <- TreePts$idTree[IndexPts]
}
}
ptsGrpl$isEmpty <- FALSE
for (h in 1:dim(ptsGrpl)[1]) {
if (ptsGrpl[h,] %>%
st_intersects(AccessMainTrails %>%
filter(ID == ID_Access),sparse = FALSE) == FALSE) {
ptsGrpl$isEmpty[h] <- TRUE
}else{
suppressWarnings(st_geometry(ptsGrpl[h,]) <- st_geometry(ptsGrpl[h,] %>%
st_intersection(AccessMainTrails %>%
filter(ID == ID_Access))))
}
}
#Filter polygons which intersect accessible area to second trails
ptsGrpl<- st_set_crs(ptsGrpl,st_crs(AccessPolygons)) # set crs from AccessPolygons
ptsGrpl <- ptsGrpl %>%
filter(isEmpty == FALSE) %>%
dplyr::select(-isEmpty) %>%
mutate(IDpts = paste0("I.", row_number()))
ptsWIP <- ptsGrpl %>% #def Grpl point as WIP point
st_set_agr("constant") %>%
st_centroid()
ptsCbl <- TreePts %>% #def cbl polygons
st_buffer(dist = advancedloggingparameters$CableLength) %>%
st_snap_to_grid(size = 1) %>%# avoid GEOS error (st-intersection issue)
st_set_precision(1) %>%
st_temporaryintersection(topography = topography,
plotmask = plotmask,
advancedloggingparameters = advancedloggingparameters) %>%
st_make_valid() %>%
mutate(type = "Inter")%>%
dplyr::select(-ID)
for (OriginI in 1:length(ptsCbl$origins)) {
for (OriginJ in 1:length(ptsCbl$origins[[OriginI]])) {
IndexPts <- ptsCbl$origins[[OriginI]][OriginJ]
ptsCbl$origins[[OriginI]][OriginJ] <- TreePts$idTree[IndexPts]
}
}
ptsCbl$isEmpty <- FALSE
for (h in 1:dim(ptsCbl)[1]) {
if (ptsCbl[h,] %>%
st_intersects(st_geometry(AccessMainTrails %>%
filter(ID == ID_Access) %>%
st_union()),sparse = FALSE) == FALSE) {
ptsCbl$isEmpty[h] <- TRUE
}else{
suppressWarnings(st_geometry(ptsCbl[h,]) <- st_geometry(ptsCbl[h,] %>%
st_intersection(st_geometry(AccessMainTrails %>%
filter(ID == ID_Access)))))
}
}
ptsCbl <- st_set_crs(ptsCbl,st_crs(AccessPolygons)) #set crs from AccessPolygons
ptsCbl <- ptsCbl %>%
filter(isEmpty == FALSE) %>%
dplyr::select(-isEmpty) %>%
mutate(IDpts = paste0("I.",row_number())) #Filter polygons which intersect accessible area to second trails
ptsWIPCbl <- ptsCbl %>%#Convert polygons to centroid
st_set_agr("constant") %>%
st_centroid()
ptsWIP <- ptsWIP %>%
arrange(desc(n.overlaps))
#Select adjacent grpl graph
SlopeCondRd <- SlopeCondGrpl
}else{
ptsCbl <- TreePts %>% #def cbl point
st_buffer(dist = advancedloggingparameters$CableLength) %>%
st_snap_to_grid(size = 1) %>% # avoid GEOS error (st-intersection issue)
st_set_precision(1) %>%
st_temporaryintersection(topography = topography,
plotmask = plotmask,
advancedloggingparameters = advancedloggingparameters) %>%
st_make_valid() %>%
mutate(IDpts = paste0("I.",row_number())) %>%
mutate(type = "Inter") %>% dplyr::select(-ID)
ptsCbl <- st_set_crs(ptsCbl,st_crs(AccessPolygons)) #set crs from AccessPolygons
for (OriginI in 1:length(ptsCbl$origins)) {
for (OriginJ in 1:length(ptsCbl$origins[[OriginI]])) {
IndexPts <- ptsCbl$origins[[OriginI]][OriginJ]
ptsCbl$origins[[OriginI]][OriginJ] <- TreePts$idTree[IndexPts]
}
}
ptsCbl$isEmpty <- FALSE
for (h in 1:dim(ptsCbl)[1]) {
if (ptsCbl[h,] %>%
st_intersects(AccessMainTrails %>%
filter(ID == ID_Access),sparse = FALSE) == FALSE) {
ptsCbl$isEmpty[h] <- TRUE
}else{
suppressWarnings(st_geometry(ptsCbl[h,]) <- st_geometry(ptsCbl[h,] %>%
st_intersection(AccessMainTrails %>%
filter(ID == ID_Access))))
}
}
ptsWIP <- ptsCbl %>%
filter(isEmpty == FALSE) %>%
dplyr::select(-isEmpty) %>%#def cbl point as WIP point
st_set_agr("constant") %>%
st_centroid()
ptsWIP <- ptsWIP %>% #filter cbl intersection centroid point out plot
filter(st_intersects(st_geometry(ptsWIP),
st_geometry(plotmask %>% st_as_sf()),
sparse = FALSE)) %>%
arrange(desc(n.overlaps))
#Select adjacent cbl graph
SlopeCondRd <- SlopeCond
}
}
RemainTree <- dim(TreePts)[1]
}
if (RemainTree > 0L){
stop("Failed to compute 2nd trails")
}
}
}
MainTrailsAccess <- AccessPointAll %>%
st_set_crs(st_crs(topography))
inventory <- treeselectionoutputs$inventory
# initial inventory
inventory0 <- inventory
lines <- do.call(rbind, Lines)
# WinchingMachine info in the inventory
TrailsIdentity_df <- as.data.frame(lines) %>%
unnest(cols = c(LineID, LoggedTrees, TypeExpl, IdMachineZone)) # unnesting flattens it back out into regular columns
TrailsIdentity_df <- TrailsIdentity_df %>%
dplyr::select(LoggedTrees, TypeExpl, IdMachineZone) %>%
rename(idTree = LoggedTrees) %>%
rename(WinchingMachine = TypeExpl) %>%
unique()
if(any(duplicated(TrailsIdentity_df$idTree))) # trees logged several times?
stop("Error: duplicated tree in logged trees")
inventory <- inventory %>%
left_join(TrailsIdentity_df, by = "idTree")
paths <- do.call(rbind, pathLines)
paths <- paths %>% st_as_sf() %>% st_make_valid() %>%
filter(!st_is_empty(paths %>% st_as_sf())) %>%
filter(st_geometry_type(.) %in% c("LINESTRING"))
if (dim(paths)[1]>0) {
paths <- paths$geometry %>% as_Spatial()
if ( length(paths)==1) {
SmoothedTrails <- paths %>% st_as_sf() %>%
st_buffer(dist = advancedloggingparameters$ScndTrailWidth/2) %>% st_make_valid() %>%
st_union()%>% st_set_crs(st_crs(topography))
TrailsDensity <- (SmoothedTrails %>%
st_intersection(plotmask %>% st_as_sf()) %>%
st_area / advancedloggingparameters$ScndTrailWidth)/(plotmask %>% st_as_sf() %>% st_area() /10000)
}else{
secondtrails <- smoothtrails(paths,
plotmask,
verbose = verbose,
advancedloggingparameters = advancedloggingparameters)
SmoothedTrails <- secondtrails$SmoothedTrails %>% st_set_crs(st_crs(topography))
TrailsDensity <- secondtrails$TrailsDensity
}
#
#
# # Records the dead trees
if(fuel == "0"){
if (!("DeathCause" %in% names(inventory))){
inventory <- inventory %>%
add_column(DeathCause = NA) # if "DeathCause" column doesnt exist create it
}
DeadTrees <- suppressWarnings(st_intersection(
st_make_valid(st_set_crs(st_as_sf(inventory, coords = c("Xutm", "Yutm")), st_crs(topography))),
st_make_valid(st_as_sf(SmoothedTrails) %>% st_union()) # "make valid" to avoid self-intersection
)) %>%
add_column(DeadTrees = "1") %>%
dplyr::select(idTree, DeadTrees) %>% st_drop_geometry() # remove geometry column (sf to data.frame)
DeadTrees <- unique(DeadTrees)
inventory <- inventory %>%
left_join(DeadTrees, by = "idTree") %>%
mutate(DeathCause = ifelse(is.na(DeathCause) & Selected != "1" & DeadTrees == "1",
"2ndtrail", DeathCause)) %>% # Damage trees
dplyr::select(-DeadTrees)
}
}else{
paths <- NULL
SmoothedTrails <- MainTrailsAccess %>% st_union() %>%
st_buffer(dist = advancedloggingparameters$MaxMainTrailWidth/2) %>% st_set_crs(st_crs(topography))
TrailsDensity <- 0
if(fuel == "0"){
if (!("DeathCause" %in% names(inventory))){
inventory <- inventory %>%
add_column(DeathCause = NA) # if "DeathCause" column doesnt exist create it
}
DeadTrees <- suppressWarnings(sf::st_intersection(
sf::st_make_valid(st_set_crs(st_as_sf(inventory, coords = c("Xutm", "Yutm")), st_crs(topography))),
sf::st_make_valid(st_as_sf(SmoothedTrails) %>% st_union()) # "make valid" to avoid self-intersection
)) %>%
add_column(DeadTrees = "1") %>%
dplyr::select(idTree, DeadTrees) %>% st_drop_geometry()
DeadTrees <- unique(DeadTrees)
inventory <- inventory %>%
left_join(DeadTrees, by = "idTree") %>%
mutate(DeathCause = ifelse(is.na(DeathCause) & Selected != "1" & DeadTrees == "1",
"2ndtrail", DeathCause)) %>% # Damage trees
dplyr::select(-DeadTrees)
}
}
if(nrow(inventory) != nrow(inventory0))
stop("The number of rows between the input inventory and the output inventory
of the function secondtrailsopening() is not the same.The function must be corrected.")
# OUTPUTS
if (WinchingInit == "2") {
secondtrails <- list("inventory" = inventory,
"SmoothedTrails" = SmoothedTrails,
"TrailsDensity" = TrailsDensity,
"TrailsIdentity" = lines, # "LineID","LoggedTrees", "TypeExpl"
"MainTrailsAccess" = MainTrailsAccess,
"RawSecondTrails" = paths,
"CostRasterAgg" = list("CostRasterMean" = CostRasterMean,"CostRasterMeanGrpl" = CostRasterMeanGrpl))
}else{
secondtrails <- list("inventory" = inventory,
"SmoothedTrails" = SmoothedTrails,
"TrailsDensity" = TrailsDensity,
"TrailsIdentity" = lines, # "LineID","LoggedTrees", "TypeExpl"
"MainTrailsAccess" = MainTrailsAccess,
"RawSecondTrails" = paths,
"CostRasterAgg" = list("CostRasterMean" = CostRasterMean,"CostRasterMeanGrpl" = NULL))
}
return(secondtrails)
}
VincyaneBadouard/LoggingLab documentation built on Oct. 16, 2024, 9:42 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions secondtrailsopening
Documented in secondtrailsopening
#' Secondary skidding trails layout
#'
#' @description Starting from the main skidding trails, draw secondary skidding
#' trails on the zones accessible to the machines, allowing to collect the
#' selected trees with the chosen machines ("winching" argument). The layout
#' is optimised to reduce the distance covered while respecting topographical
#' constraints and avoiding trees to protect.
#'
#' @param topography Digital terrain model (DTM) of the inventoried plot (LiDAR
#' or SRTM) (\code{\link{DTMParacou}}) (RasterLayer **with a crs in UTM**)
#'
#' @param plotmask Inventoried plot mask
#' (SpatialPolygonsDataFrame **with a crs in UTM**)
#'
#' @param maintrails Main trails defined at the entire harvestable area (sf
#' linestring **with a crs in UTM**)
#'
#' @param treeselectionoutputs A list with:
#' - your inventory with: "DistCriteria", "Slope", "SlopeCriteria",
#' "LoggingStatus", "Selected", "Up", "VolumeCumSum",
#' "ProbedHollowProba", "ProbedHollow" new columns
#' (see the outputs metadata in the vignette).
#'
#' - the objective volume
#'
#' - 6 sets of spatial points:
#' harvestable, selected, future and reserve, hollow and fuel wood trees
#'
#' @param harvestablepolygons Accessible zones of the inventoried plot
#' (default: \code{\link{harvestableareadefinition}}) (sfc_MULTIPOLYGON)
#'
#' @param machinepolygons Accessible zones for machines of the inventoried plot
#' (default: \code{\link{harvestableareadefinition}}) (sf polygons data.frame)
#'
#' @param maintrailsaccess Access point of maintrail for each PU (prospection
#' unit) (sf or sfc)
#'
#' @param plotslope Slopes (in radians) of the inventoried plot (with a
#' neighbourhood of 8 cells) (default:
#' \code{\link{HarvestableAreaOutputsCable}})
#' (RasterLayer **with a crs in UTM**)
#'
#'@param scenario Logging scenario among: "RIL1", "RIL2broken", "RIL2", "RIL3",
#' "RIL3fuel", "RIL3fuelhollow" or "manual"(character) (see the
#' vignette)
#'
#'@param winching
#' "0": no cable or grapple (trail to tree foot)
#' "1": only cable (default = 40m)
#' "2": grapple (default = 6m) + cable (grapple priority)
#' If grapple + cable (winching = "2") without fuel wood (fuel = "0")
#' recovery of the tree foot with grapple if possible (respected grapple
#' conditions) otherwise with cable with angle to the trail.
#' Avoidance of future/reserves if chosen.
#'
#'@param fuel Fuel wood exploitation: no exploitation = "0", exploitation of
#' damage and unused part of logged trees for fuel = "1", exploitation of
#' hollow trees, damage and and unused part of the log for fuel = "2".
#' If fuel wood harvesting is chosen, these trails are preliminary,
#' and the trails will be adjusted after tree felling for fuel wood recovery.
#' The mortality of the trails is therefore not recorded at this stage
#'
#'@param verbose Allow to provide messages from internal functions (boolean)
#'
#' @param advancedloggingparameters Other parameters of the logging simulator
#' \code{\link{loggingparameters}} (list)
#'
#' @return A list with :
#' - *inventory*: Updated inventory
#'
#' - *SmoothedTrails*: Smoothed secondary trails (MULTIPOLYGON with crs)
#'
#' - *TrailsDensity*: Secondary trails density (in m/ha)
#'
#' - *TrailsIdentity*: Information on sections of the trails (matrix) with:
#' - *LineID*:
#' - *LoggedTrees*: idTree of trees reached by the trails
#' - *TypeExpl*: type of winching
#'
#' - *MainTrailsAccess*: Random access point of main trail for each harvestable
#' zone (sfc_POINT with crs)
#'
#' - *RawSecondTrails*: Non-smoothed secondary trails (SpatialLines with crs)
#'
#' - *CostRasterAgg*: The cost raster (RasterLayer with crs)
#'
#' @importFrom sf st_cast st_as_sf st_as_sfc st_intersection st_union st_sample st_join
#' st_buffer as_Spatial st_centroid st_set_precision st_make_valid st_set_agr
#' st_geometry st_area st_is_empty st_set_crs st_crs sf_use_s2 st_geometry<-
#' st_drop_geometry st_geometry_type
#'
#' @importFrom dplyr mutate row_number select as_tibble left_join if_else filter
#' arrange desc distinct
#'
#' @importFrom raster raster extend extent focal res crs mask crop rasterize
#' rasterToPoints rasterToPolygons rasterFromXYZ aggregate values ncell
#' values<- adjacent resample ncol
#'
#' @importFrom dplyr select filter mutate
#' @importFrom tidyr unnest
#' @importFrom tibble add_column
#' @importFrom sp proj4string<- coordinates<-
#' @importFrom sf st_as_sf st_union st_distance st_crs
#' @importFrom lwgeom st_snap_to_grid
#' @importFrom smoothr smooth
#' @importFrom gdistance transition geoCorrection
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @importFrom stats na.exclude
#'
#' @export
#'
#' @examples
#' set.seed(1)
#' data(Paracou6_2016)
#' data(DTMParacou)
#' data(PlotMask)
#' data(SpeciesCriteria)
#' data(HarvestableAreaOutputsCable)
#' data(MainTrails)
#'
#' scenario <- "manual"
#' winching <- "2"
#' objective <- 10
#' fuel <- "2"
#' diversification <- TRUE
#'
#' inventory <- addtreedim(cleaninventory(Paracou6_2016, PlotMask),
#' volumeparameters = ForestZoneVolumeParametersTable)
#'
#' treeselectionoutputs <- suppressWarnings(treeselection(inventory,
#' topography = DTMParacou,
#' speciescriteria = SpeciesCriteria,
#' scenario = "manual", objective = objective,
#' fuel = fuel,
#' diversification = diversification,
#' winching = winching,
#' specieslax = FALSE, objectivelax = TRUE,
#' harvestablearea = HarvestableAreaOutputsCable$HarvestableArea,
#' plotslope = HarvestableAreaOutputsCable$PlotSlope,
#' maintrails = MainTrails,
#' harvestablepolygons = HarvestableAreaOutputsCable$HarvestablePolygons,
#' advancedloggingparameters = loggingparameters()))
#'
#' \dontrun{
#' secondtrails <- secondtrailsopening(
#' topography = DTMParacou,
#' plotmask = PlotMask,
#' maintrails = MainTrails,
#' plotslope = HarvestableAreaOutputsCable$PlotSlope,
#' harvestablepolygons = HarvestableAreaOutputsCable$HarvestablePolygons,
#' machinepolygons = HarvestableAreaOutputsCable$MachinePolygons,
#' treeselectionoutputs = treeselectionoutputs,
#' scenario = "manual",
#' winching = winching,
#' fuel = "0",
#' advancedloggingparameters = loggingparameters())
#' }
#'
#' data(SecondaryTrails)
#' library(ggplot2)
#' library(sf)
#'
#' NewInventory <- SecondaryTrails$inventory
#'
#'
#' Harvestable <- sf::st_as_sf(
#' dplyr::filter(NewInventory, LoggingStatus == "harvestable"),
#' coords = c("Xutm", "Yutm")) %>%
#' st_set_crs(st_crs(SecondaryTrails$MainTrailsAccess))
#'
#' HarvestableUp <- sf::st_as_sf(
#' dplyr::filter(NewInventory, LoggingStatus == "harvestableUp"),
#' coords = c("Xutm", "Yutm")) %>%
#' st_set_crs(st_crs(SecondaryTrails$MainTrailsAccess))
#'
#' Selected <- sf::st_as_sf(
#' dplyr::filter(NewInventory, Selected == "1"), coords = c("Xutm", "Yutm")) %>%
#' st_set_crs(st_crs(SecondaryTrails$MainTrailsAccess))
#'
#' Reserve <- sf::st_as_sf(
#' dplyr::filter(NewInventory, LoggingStatus == "reserve"),
#' coords = c("Xutm", "Yutm")) %>%
#' st_set_crs(st_crs(SecondaryTrails$MainTrailsAccess))
#'
#' Future <- sf::st_as_sf(
#' dplyr::filter(NewInventory, LoggingStatus == "future"),
#' coords = c("Xutm", "Yutm")) %>%
#' st_set_crs(st_crs(SecondaryTrails$MainTrailsAccess))
#'
#' ProbedHollow <- sf::st_as_sf(
#' dplyr::filter(NewInventory, ProbedHollow == "1"),
#' coords = c("Xutm", "Yutm")) %>%
#' st_set_crs(st_crs(SecondaryTrails$MainTrailsAccess))
#'
#' VisibleDefect <- sf::st_as_sf(
#' dplyr::filter(NewInventory, VisibleDefect == "1"),
#' coords = c("Xutm", "Yutm")) %>%
#' st_set_crs(st_crs(SecondaryTrails$MainTrailsAccess))
#'
#'
#' ggplot() +
#' # Harvestable zones
#' geom_sf(data = HarvestableAreaOutputsCable$HarvestablePolygons,
#' fill = "olivedrab", alpha = 0.1) +
#' geom_sf(data = HarvestableAreaOutputsCable$MachinePolygons,
#' fill = "olivedrab", alpha = 0.5) +
#' labs(alpha = "Harvestable") +
#' labs(title = "Paracou P6 - secondary trails") +
#'
#' geom_sf(data = VisibleDefect,
#' aes(colour = "Visible defect"), show.legend = "point") +
#' # geom_sf(data = Future,
#' # aes(colour = "Future"), show.legend = "point", size = 4) +
#' geom_sf(data = Reserve,
#' aes(colour = "Reserve"), show.legend = "point", size = 4) +
#' geom_sf(data = Harvestable,
#' aes(colour = "Harvestable"), show.legend = "point", size = 4) +
#' geom_sf(data = HarvestableUp,
#' aes(colour = "HarvestableUp"), show.legend = "point", size = 4) +
#' geom_sf(data = Selected,
#' aes(colour = "Selected"), show.legend = "point") +
#' geom_sf(data = ProbedHollow,
#' aes(colour = "Probed hollow"), show.legend = "point") +
#'
#' # 2ndary trails
#' geom_sf(data = st_as_sf(SecondaryTrails$SmoothedTrails), col = "darkgreen") +
#' geom_sf(data = st_as_sf(SecondaryTrails$MainTrailsAccess), col = "black") +
#'
#' scale_colour_manual(values = c("Visible defect" = "pink",
#' "Harvestable" = "skyblue", "HarvestableUp" = "blue", "Selected" = "red",
#' "Future" = "orange", "Reserve" = "purple", "Probed hollow" = "forestgreen",
#' "Second trails" = "darkgreen", "Harvestable area" = "olivedrab"))
#'
#' SecondaryTrails$TrailsIdentity
#'
secondtrailsopening <- function(
topography,
plotmask,
maintrails,
plotslope,
harvestablepolygons,
machinepolygons,
maintrailsaccess = NULL,
treeselectionoutputs,
scenario,
winching = NULL,
fuel = NULL,
verbose = FALSE,
advancedloggingparameters = loggingparameters()
){
# Arguments check
if(!inherits(treeselectionoutputs, "list"))
stop("The 'treeselectionoutputs' arguments of the 'secondtrailsopening'
function must be list following treeselection output format")
if(!inherits(plotmask, "SpatialPolygonsDataFrame"))
stop("The 'plotmask' argument of the 'secondtrailsopening' function must be
SpatialPolygonsDataFrame")
# if(!any(unlist(lapply(list(maintrails), inherits, "sf" ))))
# stop("The 'maintrails' argument of the 'secondtrailsopening' function must be sf polygon")
if(!inherits(topography, "RasterLayer"))
stop("The 'topography' argument of the 'secondtrailsopening' function must
be RasterLayer")
if(sf::st_is_empty(treeselectionoutputs$SelectedTreesPoints)[1]) {
stop("The 'treeselectionoutputs' argument does not contain any selected
tree.")
}
# Options
options("rgdal_show_exportToProj4_warnings"="none") # to avoid gdal warnings
# gc() # remove intermediary files (to avoid gdal warnings)
# Global Variables
slope <- x <- y <- Harvestable <- idTree <- ID <- type <- ptAcc <- . <- NULL
EstCost <- n.overlaps <- TypeAcc <- IDpts <- Logged <- AccessPolygons <- NULL
Selected <- DeathCause <- ID_Acc <- isEmpty <- gprlAcc <- cblAcc <- NULL
ID.y <- IdMachineZone <- IdMachineZone.y <- IdMachineZone.x <- LineID <- LoggedTrees <- TypeExpl <- NULL
#### Redefinition of the parameters according to the chosen scenario ####
scenariosparameters <- scenariosparameters(scenario = scenario, winching = winching, fuel = fuel)
WinchingInit <- scenariosparameters$winching
winching <- WinchingInit
fuel <- scenariosparameters$fuel
##################################
suppressMessages(sf_use_s2(FALSE)) # to deal with actual unresolved s2 issues in sf ("Spherical geometry (s2) switched off")
CostMatrix <- advancedloggingparameters$CostMatrix
factagg <- floor(advancedloggingparameters$SlopeDistance/res(topography)[1])
advancedloggingparameters$CableLength <- advancedloggingparameters$CableLength + factagg
# Transformation of the DTM so that the maintrails are outside the plot
DTMExtExtended <- raster::extend(topography, c(2*factagg,2*factagg)) # extend the extent
fill.boundaries <- function(x) {
center = 0.5 + (factagg*factagg/2)
if( is.na(x)[center] ) {
return( round(mean(x, na.rm=T),5) )
} else {
return( x[center] )
}
}
DTMExtended <- raster::focal(DTMExtExtended,
matrix(1,factagg,factagg),
fun=fill.boundaries,
na.rm=F, pad=T)
for (iterExt in 1:5) {
DTMExtended <- raster::focal(DTMExtended,
matrix(1,factagg,factagg),
fun=fill.boundaries,
na.rm=F, pad=T)
}
# Generate accessible area from HarvestablePolygones and winching choice
AccessPolygons <- machinepolygons
maintrailsRed <- maintrails %>% st_as_sfc() %>%
st_cast("POLYGON") %>% st_buffer(-factagg) %>%
st_cast("LINESTRING") %>% st_buffer(factagg+0.5)
#-------------------------------------------
if (!is.null(maintrailsaccess)) {
AccessMainTrails <- AccessPolygons %>% st_union() %>%
st_cast("POLYGON", warn = FALSE) %>%
st_as_sf() %>% st_join(maintrailsaccess) %>%
select(ID)
PartMainTrails <- st_intersection(st_geometry(maintrailsRed),
st_geometry(AccessMainTrails)) %>%
st_union(by_feature = T) %>% st_as_sf() %>%
filter(st_geometry_type(.) %in% c("POLYGON","MULTIPOLYGON") ) %>%
st_cast("MULTIPOLYGON", warn = FALSE) %>% # "In st_cast.MULTIPOLYGON(X[[i]], ...) : polygon from first part only"
st_as_sf() %>%
st_set_agr(value = "constant") %>%
st_join(AccessMainTrails)%>%
st_make_valid()
PartMainTrails <- PartMainTrails %>%
arrange(desc(st_area(PartMainTrails))) %>%
filter(duplicated(PartMainTrails$ID) == FALSE)
AccessPointAll <- maintrailsaccess
}else{
accesspts <- genaccesspts(topography = topography,
machinepolygons = machinepolygons,
maintrails = maintrails,
advancedloggingparameters = advancedloggingparameters)
PartMainTrails <- accesspts$PartMainTrails
# Generate point access in the intersections between accessible area and maintrails (ID = accessible area index)
AccessPointAll <- accesspts$AccessPointAll
AccessMainTrails <- AccessPolygons %>% st_union() %>%
st_cast("POLYGON", warn = FALSE) %>%
st_as_sf() %>% st_join(accesspts$AccessPointAll) %>%
select(ID)
}
#-------------------------------------------
# Generate Cost raster --> cf CostMatrix
CostRaster <- raster(extent(DTMExtended),resolution = res(DTMExtended), crs = crs(DTMExtended))
values(CostRaster) <- CostMatrix[[2]][[1]]$CostValue
#Generate weight according to slope conditions
# RastersToPoints
plotslopePoint <- as_tibble(rasterToPoints(plotslope))
CostRasterPoint <- as_tibble(rasterToPoints(CostRaster))
# left join par x et y
plotTib <-
left_join(plotslopePoint, CostRasterPoint, by = c('x', 'y'))
CostSlopeRaster <- plotTib %>%
mutate(Harvestable = if_else(
slope <= atan(CostMatrix[[1]][[1]]$Slope/100),
true = CostMatrix[[1]][[1]]$Cost,
false = if_else(
slope <= atan(CostMatrix[[1]][[2]]$Slope/100),
true = CostMatrix[[1]][[2]]$Cost,
false = if_else(
slope <= atan(CostMatrix[[1]][[3]]$Slope/100),
true = CostMatrix[[1]][[3]]$Cost,
false = if_else(
slope <= atan(CostMatrix[[1]][[4]]$Slope/100),
true = CostMatrix[[1]][[4]]$Cost,
false = if_else(
slope <= atan(CostMatrix[[1]][[5]]$Slope/100),
true = CostMatrix[[1]][[5]]$Cost,
false = CostMatrix[[1]][[6]]$Cost)
)
)
)
)) %>% dplyr::select(x,y,Harvestable)
CostRaster <- rasterFromXYZ(CostSlopeRaster, crs = crs(DTMExtended)) # Update weights from plotTib tibble
#Aggregation each raster to selected resolution
CostRasterMean <- aggregate(CostRaster, fact=factagg, fun=mean)
DTMmean <- aggregate(DTMExtended, fact=factagg, fun=mean)
plotmaskagg <- rasterToPolygons(DTMmean, dissolve=T) %>% st_as_sf() %>% st_union() %>% as_Spatial()
#Generate protection buffer for Big Trees (dist : ScndTrailWidth/2 + 2 m)
BigTreesRaster <- raster(extent(DTMmean),resolution = res(DTMmean), crs = crs(DTMmean))
values(BigTreesRaster) <- 0
BigTreesRaster <- rasterize(x = as_Spatial(treeselectionoutputs$BigTreesPoints %>%
st_buffer(dist = advancedloggingparameters$ScndTrailWidth/2 + 2) ),
y = BigTreesRaster,
field = CostMatrix[[2]][[3]]$CostValue,
update = TRUE)
BigTreesRaster <- crop(BigTreesRaster, CostRasterMean)
BigTreesRaster <- mask(BigTreesRaster, plotmaskagg)
#Generate protection buffer for Reserve Trees (dist : ScndTrailWidth/2 + 2 m)
ReserveRaster <- raster(extent(DTMmean),resolution = res(DTMmean), crs = crs(DTMmean))
values(ReserveRaster) <- 0
ReserveRaster <- rasterize(x = as_Spatial(treeselectionoutputs$HarvestableTreesPoints %>%
st_buffer(dist = advancedloggingparameters$ScndTrailWidth/2 + 2) ),
y = ReserveRaster,
field = CostMatrix[[2]][[4]]$CostValue,
update = TRUE)
ReserveRaster <- crop(ReserveRaster, CostRasterMean)
ReserveRaster <- mask(ReserveRaster, plotmaskagg)
#Generate protection buffer for Futures Trees (dist : ScndTrailWidth/2 + 2 m)
FutureRaster <- raster(extent(DTMmean),resolution = res(DTMmean), crs = crs(DTMmean))
values(FutureRaster ) <- 0
FutureRaster <- rasterize(x = as_Spatial(treeselectionoutputs$HarvestableTreesPoints %>%
st_buffer(dist = advancedloggingparameters$ScndTrailWidth/2 +2) ),
y = FutureRaster ,
field = CostMatrix[[2]][[5]]$CostValue,
update = TRUE)
FutureRaster <- crop(FutureRaster, CostRasterMean)
FutureRaster <- mask(FutureRaster, plotmaskagg)
#Update Cost Raster with protection buffers
#Generate protection buffer for selected Trees (dist : ScndTrailWidth/2 + 2 m)
SelectedRaster <- raster(extent(DTMmean),resolution = res(DTMmean), crs = crs(DTMmean))
values(SelectedRaster ) <- 0
SelectedRaster <- rasterize(x = as_Spatial(treeselectionoutputs$SelectedTreesPoints %>%
st_buffer(dist = 2) ),
y = SelectedRaster,
field = CostMatrix[[2]][[3]]$CostValue/2,
update = TRUE)
SelectedRaster <- crop(SelectedRaster, CostRasterMean)
SelectedRaster <- mask(SelectedRaster, plotmaskagg)
#Update Cost Raster with protection buffers
CostRasterMean <- CostRasterMean + BigTreesRaster + ReserveRaster + FutureRaster + SelectedRaster
#Generate maintrail intersection cost raster
CostRasterMean <- rasterize(x = as_Spatial(AccessPointAll %>% st_buffer(dist = advancedloggingparameters$ScndTrailWidth+2)),
y = CostRasterMean ,
field = CostMatrix[[2]][[6]]$CostValue,
update = TRUE)
#Generate Slope accessibility for grapple machine
if (winching == "2") {
plotslopePointGrpl <- as_tibble(rasterToPoints(plotslope))
CostRasterPointGrpl <- as_tibble(rasterToPoints(CostRaster))
# left join par x et y
plotTibGrpl <-
left_join(plotslopePointGrpl, CostRasterPointGrpl, by = c('x', 'y'))
CostSlopeRasterGrpl <- plotTibGrpl %>%
mutate(Harvestable = if_else(
slope <= atan(advancedloggingparameters$GrappleMaxslope/100),
true = 0,
false = Inf
)) %>% dplyr::select(x,y,Harvestable)
CostRasterGrpl <- rasterFromXYZ(CostSlopeRasterGrpl, crs = crs(CostRasterMean))
PolygonGrpl <- rasterToPolygons(x = CostRasterGrpl,
n = 8,
dissolve = TRUE)
CostRasterMeanGrpl <- rasterize(x = PolygonGrpl[PolygonGrpl$Harvestable == Inf,],
y = CostRasterMean ,
field = Inf,
update = TRUE)
CostRasterMeanGrpl <- rasterize(x = as_Spatial(AccessPointAll %>% st_buffer(dist = advancedloggingparameters$ScndTrailWidth+2)),
y = CostRasterMeanGrpl ,
field = CostMatrix[[2]][[6]]$CostValue,
update = TRUE)
#Generate accessible weights raster
AccessRaster <- raster(extent(CostRaster),resolution = res(CostRaster), crs = crs(CostRaster))
values(AccessRaster) <- CostMatrix[[2]][[2]]$CostValue
AccessRaster <- rasterize(x = as_Spatial(machinepolygons %>% st_buffer(dist = factagg)),
y = AccessRaster ,
field = 1,
update = TRUE)
AccessRaster <- aggregate(AccessRaster,fact=factagg, fun=min)
AccessRaster <- crop(AccessRaster, CostRasterMeanGrpl)
AccessRaster <- mask(AccessRaster, plotmaskagg)
#Update Cost Raster with accessible weights raster
CostRasterMeanGrpl <- CostRasterMeanGrpl * AccessRaster
}
#Generate accessible weights raster
AccessRaster <- raster(extent(CostRaster),resolution = res(CostRaster), crs = crs(CostRaster))
values(AccessRaster) <- CostMatrix[[2]][[2]]$CostValue
AccessRaster <- rasterize(x = as_Spatial(machinepolygons %>% st_buffer(dist = factagg)),
y = AccessRaster ,
field = 1,
update = TRUE)
AccessRaster <- aggregate(AccessRaster,fact=factagg, fun=min)
AccessRaster <- crop(AccessRaster, CostRasterMean)
AccessRaster <- mask(AccessRaster,plotmaskagg)
#Update Cost Raster with accessible weights raster
CostRasterMean <- CostRasterMean * AccessRaster
pathLines <- list() #Initialize storage pathlines
Lines <- list() #Initialize storage logged trees
k <- 1 #Initialize pathlines index
l <- 1 #Initialize lines index
#Generate appropriate selected trees points format
ptsAllinit <- treeselectionoutputs$SelectedTreesPoints %>%
dplyr::select(idTree) %>%
st_cast("POINT", warn = FALSE) %>%
mutate(ID = NA) %>%
mutate(type = "Tree") %>%
st_set_crs(st_crs(AccessPointAll)) %>%
st_join(AccessMainTrails %>% st_buffer(dist = max((advancedloggingparameters$ScndTrailWidth/2),factagg)) %>%
mutate(ID_Acc = ID) %>% dplyr::select(-ID)) %>%
dplyr::select(ID,ID_Acc,type,idTree)
ptsAllNA <- ptsAllinit %>%
filter(is.na(ID_Acc)) %>%
mutate(ID_Acc = AccessMainTrails$ID[c(max.col(-st_distance(ptsAllinit %>%
filter(is.na(ID_Acc)),AccessMainTrails)))])
ptsAllNotNA <- ptsAllinit %>%
filter(!is.na(ID_Acc))
ptsAll <- rbind(ptsAllNotNA,ptsAllNA)
# Reassign Selected Tree values (= BigTrees) to the aggregated Cost raster
CostRasterMean <- rasterize(x = as_Spatial(ptsAll %>% st_buffer(dist = factagg)),
y = CostRasterMean ,
field = CostMatrix[[2]][[3]]$CostValue,
update = TRUE)
if (winching == "2") {
# Reassign Selected Tree values (= BigTrees) to the aggregated Cost raster (Grpl)
CostRasterMeanGrpl <- rasterize(x = as_Spatial(ptsAll %>%
st_buffer(dist = factagg)),
y = CostRasterMeanGrpl ,
field = CostMatrix[[2]][[3]]$CostValue,
update = TRUE)
CostRasterMeanGrpl <- crop(CostRasterMeanGrpl, CostRasterMean)
CostRasterMeanGrpl <- mask(CostRasterMeanGrpl, plotmaskagg)
}
#Generate maintrail intersection cost raster
CostRasterMean <- rasterize(x = as_Spatial(AccessPointAll %>% st_buffer(dist = 3*factagg)),
y = CostRasterMean ,
field = CostMatrix[[2]][[6]]$CostValue,
update = TRUE)
CostRasterMean <- crop(CostRasterMean,DTMmean)
CostRasterMean <- mask(CostRasterMean,plotmaskagg)
DTMmean <- crop(DTMmean,CostRasterMean)
DTMmean <- mask(DTMmean,plotmaskagg)
#Compute conductance raster
CondSurf <- 1/CostRasterMean
########### Compute slope criteria transition graph ###############
if (winching == "2") {
#Compute adjacent transition layer according to slope conditions (winching = "2")
SlopeCondGrpl <- sloperdcond(topography = DTMmean,
advancedloggingparameters = advancedloggingparameters,
grapple = TRUE)
}
#Compute adjacent transition layer according to slope conditions (winching = "1")
SlopeCond <- sloperdcond(topography = DTMmean,advancedloggingparameters = advancedloggingparameters)
########### Compute LCP algorithm ###############
for(ID_Access in unique(ptsAll$ID_Acc)) { # ici
winching <- WinchingInit
pts <- ptsAll %>%
filter(ID_Acc == ID_Access) %>% dplyr::mutate(ID = ID_Access) %>%
dplyr::select(-ID_Acc) %>%
dplyr::distinct()
AccessPoint <- AccessPointAll %>%
filter(ID == ID_Access)
pathLinesWIP <- c()
ki <- 1
if (winching == "0") {
# winching 0 #########
#Update Cost Raster with accessible weights raster
CostRasterMean <- CostRasterMean * AccessRaster
CondSurf <- 1/CostRasterMean
pts <- st_set_crs(pts, st_crs(AccessPoint))# set crs
pts <- rbind(AccessPoint %>%
mutate(type = "Access") %>%
mutate(IPpts = paste0("A.",row_number())),
pts %>%
mutate(IPpts = paste0("T.",idTree)))
PointAcc <- pts %>%
filter(type == "Access") %>%
mutate(EstCost = NA)
#Compute Cost between points and Access points
CostDistEst <- adjtopolcp(costSurface = CondSurf,
topography = DTMmean ,
pts = as_Spatial(pts),
slopeRdCond = SlopeCond,
paths = FALSE) [,1:dim(AccessPoint)[1]]
CostDistEst <- CostDistEst[(length(PointAcc$ID)+1):length(CostDistEst)[1]]
CostDistEst <- matrix(CostDistEst,ncol = length(PointAcc$ID) )
#Attribute a least cost point access to each points
PointTree <- pts %>% filter(type == "Tree") %>%
mutate(ptAcc = max.col(-CostDistEst))
PointTree$ID <- as.factor(PointTree$ptAcc)
levels(PointTree$ID) <- as.character(AccessPoint$ID)
pts <- rbind(AccessPoint %>%
mutate(type = "Access") %>%
mutate(ptAcc = row_number()) %>%
mutate(IPpts = paste0("A.",row_number())),
PointTree %>%
mutate(IPpts = paste0("T.",idTree)))
PointAcc <- pts %>%
filter(type == "Access") %>%
mutate(EstCost = NA)
#Attribute a least cost point access to each points
PointTree <- pts %>% filter(type == "Tree") %>%
mutate(ptAcc = max.col(-CostDistEst))
if (length(PointAcc$ID) == 1) {
PointTree <- PointTree %>% mutate(EstCost = CostDistEst)
}else{PointTree <- PointTree %>% mutate(EstCost = apply(CostDistEst,1, min))}
selectedPoints <- rbind(PointAcc,PointTree)
accessPtId <- 1
TmpSelectedPts <- selectedPoints %>%
filter(ptAcc == accessPtId)
TmpAccPts <- TmpSelectedPts %>%
filter(type == "Access")
#Select Costliest selected tree
TmpTreePts <- TmpSelectedPts %>%
filter(type == "Tree") %>%
arrange(desc(EstCost))
if (dim(TmpTreePts)[1] != 0) {
for (TreeId in 1:dim(TmpTreePts)[1]){
#Compute Least cost path polygons to the WIP selected tree
TmpPtsWIP <- rbind(TmpAccPts,TmpTreePts[TreeId,])
TmpPathWIP <- adjtopolcp(costSurface = CondSurf,topography = DTMmean , pts = TmpPtsWIP %>% as_Spatial(),
slopeRdCond = SlopeCond,paths = TRUE)
if (TmpPathWIP[[1]][1,2] == 0) {
#Update Cost raster with LCP second trail
CostRasterMean <- rasterize(x = TmpTreePts[TreeId,] ,
y = CostRasterMean ,
field = CostMatrix[[2]][[7]]$CostValue,update = TRUE)
}else{
#Update Cost raster with LCP second trail
CostRasterMean <- rasterize(x = TmpPathWIP[[2]] ,
y = CostRasterMean ,
field = CostMatrix[[2]][[7]]$CostValue,update = TRUE)
}
CondSurf <- 1/CostRasterMean
#Store pathline
pathLines[[k]] <- TmpPathWIP[[2]]
pathLines[[k]]@lines[[1]]@ID <- paste("Path", TmpPtsWIP$idTree[2], sep = ".")
Lines[[k]] <- list("LineID" = k,"LoggedTrees" = TmpPtsWIP$idTree[2],"TypeExpl" = "FoT")
k <- k +1
}
}
}else{
# winching 1/2 #########
# Select intersection points from buffer polygons
PointAcc <- AccessPoint %>% #def Access Point
mutate(type = "Access") %>%
mutate(IDpts = paste0("A.",row_number())) %>%
mutate(n.overlaps = NA, origins = NA) %>%
dplyr::select(-ID)
TreePts <- pts %>% filter(type == "Tree")
if (winching == "2") {
ptsGrpl <- TreePts %>% #def Grpl point
st_buffer(dist = advancedloggingparameters$GrappleLength) %>%
st_snap_to_grid(size = .2) %>% # avoid GEOS error (st-intersection issue)
#st_set_precision(1) %>%
st_temporaryintersection(topography = topography,
plotmask = plotmask,
advancedloggingparameters = advancedloggingparameters) %>%
st_make_valid() %>%
mutate(type = "Inter") %>%
dplyr::select(-ID)
for (OriginI in 1:length(ptsGrpl$origins)) {
for (OriginJ in 1:length(ptsGrpl$origins[[OriginI]])) {
IndexPts <- ptsGrpl$origins[[OriginI]][OriginJ]
ptsGrpl$origins[[OriginI]][OriginJ] <- TreePts$idTree[IndexPts]
}
}
ptsGrpl$isEmpty <- FALSE
for (h in 1:dim(ptsGrpl)[1]) {
if (ptsGrpl[h,] %>% st_intersects(AccessMainTrails %>%
filter(ID == ID_Access),sparse = FALSE) == FALSE) {
ptsGrpl$isEmpty[h] <- TRUE
}else{
suppressWarnings(st_geometry(ptsGrpl[h,]) <- st_geometry(ptsGrpl[h,] %>% st_intersection(AccessMainTrails %>%
filter(ID == ID_Access))))
}
}
#Filter polygons which intersect accessible area to second trails
ptsGrpl<- st_set_crs(ptsGrpl,st_crs(AccessPolygons)) # set crs from AccessPolygons
ptsGrpl <- ptsGrpl %>%
filter(isEmpty == FALSE) %>%
dplyr::select(-isEmpty) %>%
mutate(IDpts = paste0("I.", row_number()))
ptsWIP <- ptsGrpl %>% #def Grpl point as WIP point
st_set_agr("constant") %>%
st_centroid()
ptsCbl <- TreePts %>% #def cbl polygons
st_buffer(dist = advancedloggingparameters$CableLength) %>%
st_snap_to_grid(size = 1) %>%# avoid GEOS error (st-intersection issue)
st_set_precision(1) %>%
st_temporaryintersection(topography = topography,
plotmask = plotmask,
advancedloggingparameters = advancedloggingparameters) %>%
st_make_valid() %>%
mutate(type = "Inter")%>%
dplyr::select(-ID)
for (OriginI in 1:length(ptsCbl$origins)) {
for (OriginJ in 1:length(ptsCbl$origins[[OriginI]])) {
IndexPts <- ptsCbl$origins[[OriginI]][OriginJ]
ptsCbl$origins[[OriginI]][OriginJ] <- TreePts$idTree[IndexPts]
}
}
ptsCbl$isEmpty <- FALSE
for (h in 1:dim(ptsCbl)[1]) {
if (ptsCbl[h,] %>% st_intersects(AccessMainTrails %>%
filter(ID == ID_Access),sparse = FALSE) == FALSE) {
ptsCbl$isEmpty[h] <- TRUE
}else{
suppressWarnings(st_geometry(ptsCbl[h,]) <- st_geometry(ptsCbl[h,] %>%
st_intersection(AccessMainTrails %>%
filter(ID == ID_Access))))
}
}
ptsCbl <- st_set_crs(ptsCbl,st_crs(AccessPolygons)) #set crs from AccessPolygons
ptsCbl <- ptsCbl %>% filter(isEmpty == FALSE) %>% dplyr::select(-isEmpty) %>%
mutate(IDpts = paste0("I.",row_number())) #Filter polygons which intersect accessible area to second trails
ptsWIPCbl <- ptsCbl %>%#Convert polygons to centroid
st_set_agr("constant") %>%
st_centroid()
ptsWIP <- ptsWIP %>%
arrange(desc(n.overlaps))
RemainTree <- max(length(unique(unlist(ptsWIP$origins))),length(unique(unlist(ptsWIPCbl$origins))))
#Select adjacent grpl graph
SlopeCondRd <- SlopeCondGrpl
}else{
ptsCbl <- TreePts %>% #def cbl point
st_buffer(dist = advancedloggingparameters$CableLength) %>%
st_snap_to_grid(size = 1) %>% # avoid GEOS error (st-intersection issue)
st_set_precision(1) %>%
st_temporaryintersection(topography = topography,
plotmask = plotmask,
advancedloggingparameters = advancedloggingparameters) %>%
st_make_valid() %>%
mutate(IDpts = paste0("I.",row_number())) %>%
mutate(type = "Inter") %>% dplyr::select(-ID)
ptsCbl <- st_set_crs(ptsCbl,st_crs(AccessPolygons)) #set crs from AccessPolygons
for (OriginI in 1:length(ptsCbl$origins)) {
for (OriginJ in 1:length(ptsCbl$origins[[OriginI]])) {
IndexPts <- ptsCbl$origins[[OriginI]][OriginJ]
ptsCbl$origins[[OriginI]][OriginJ] <- TreePts$idTree[IndexPts]
}
}
ptsCbl$isEmpty <- FALSE
for (h in 1:dim(ptsCbl)[1]) {
if (st_geometry(ptsCbl[h,]) %>% st_intersects(AccessMainTrails %>%
filter(ID == ID_Access),sparse = FALSE) == FALSE) {
ptsCbl$isEmpty[h] <- TRUE
}else{
suppressWarnings(st_geometry(ptsCbl[h,]) <- st_geometry(ptsCbl[h,] %>%
st_intersection(AccessMainTrails %>%
filter(ID == ID_Access))))
ptsCbl$isEmpty[h] <- FALSE
}
}
ptsWIP <- ptsCbl %>% filter(isEmpty == FALSE) %>%
dplyr::select(-isEmpty) %>% #def cbl point as WIP point
st_set_agr("constant") %>%
st_centroid()
ptsWIP <- ptsWIP %>% #filter cbl intersection centroid point out plot
filter(st_intersects(st_geometry(ptsWIP),
st_geometry(plotmask %>% st_as_sf()),
sparse = FALSE)) %>%
arrange(desc(n.overlaps))
RemainTree <- length(unique(ptsWIP$idTree))
#Select adjacent cbl graph
SlopeCondRd <- SlopeCond
}
#Define a switch from grapple accessible tree exploitation to cable exploitation
if (winching == "2") {
WinchingSwitch <- FALSE
}else{
WinchingSwitch <- TRUE
}
steps <- 0
RemainTreeInit <- RemainTree
#Loop over possible intersection
while (RemainTree !=0L & steps < 5*RemainTreeInit) {
steps <- steps +1
if (verbose) {
print(RemainTree)
}
#Switch from grpl to cbl exploitation when grapple accessible tree != 0
Grpl2CblFlag <- FALSE
if (length(unique(unlist(ptsWIP$origins))) == 0 & WinchingSwitch == FALSE & winching == "2") {
winching <- "1"
ptsWIP <- ptsWIPCbl
WinchingSwitch <- TRUE
}
ptsWIPmax <- rbind(PointAcc,ptsWIP %>%
filter(st_intersects(.,st_geometry(plotmask %>%
st_as_sf()),sparse = FALSE)) %>%
filter(n.overlaps == max(ptsWIP$n.overlaps))) %>%
mutate(TypeAcc = NA) %>%
mutate(EstCost = NA) %>%
mutate(ptsAcc = NA)
ptsTreeWIP <- rbind(PointAcc,TreePts %>%
mutate(n.overlaps = NA, origins = idTree,IDpts = NA) %>%
dplyr::select(-ID)) %>%
mutate(TypeAcc = NA) %>%
mutate(EstCost = NA) %>%
mutate(ptsAcc = NA)
if (ki == 1 ) {
if (winching == "2") {
ptsDirAcc <- ptsTreeWIP %>%
mutate(gprlAcc = c(FALSE,as.numeric(st_distance(ptsTreeWIP)[2:dim(ptsTreeWIP)[1],1]) < advancedloggingparameters$GrappleLength)) %>%
filter(gprlAcc == TRUE | type == "Access") %>%
dplyr::select(-gprlAcc)
TmpTypeAcc <- "Grpl"
}else{
ptsDirAcc <- ptsTreeWIP %>%
mutate(cblAcc = c(FALSE,as.numeric(st_distance(ptsTreeWIP)[2:dim(ptsTreeWIP)[1],1]) < advancedloggingparameters$CableLength)) %>%
filter(cblAcc == TRUE | type == "Access") %>%
dplyr::select(-cblAcc)
TmpTypeAcc <- "Cbl"
}
}else{
# Set pathsWIP CRS
pathsWIP <- st_set_crs(pathsWIP, st_crs(PointAcc))
if (winching == "2") {
ptsDirAcc <- ptsTreeWIP %>%
mutate(gprlAcc = c(FALSE,as.numeric(st_distance(ptsTreeWIP,st_union(pathsWIP %>%
st_buffer(dist = advancedloggingparameters$MaxMainTrailWidth/2) %>%
st_make_valid() %>%
st_union(),PointAcc %>%
st_buffer(dist = advancedloggingparameters$MaxMainTrailWidth/2) %>%
st_make_valid() %>%
st_union()))[2:dim(ptsTreeWIP)[1],1]) < advancedloggingparameters$GrappleLength)) %>%
filter(gprlAcc == TRUE | type == "Access") %>%
dplyr::select(-gprlAcc)
TmpTypeAcc <- "Grpl"
}else{
ptsDirAcc <- ptsTreeWIP %>% mutate(cblAcc = c(FALSE,as.numeric(st_distance(ptsTreeWIP,st_union(pathsWIP %>% st_buffer(dist = advancedloggingparameters$MaxMainTrailWidth/2) %>% st_make_valid() %>%
st_union(),PointAcc %>% st_buffer(dist = advancedloggingparameters$MaxMainTrailWidth/2) %>% st_make_valid() %>%
st_union()))[2:dim(ptsTreeWIP)[1],1]) < advancedloggingparameters$CableLength)) %>%
filter(cblAcc == TRUE | type == "Access") %>% dplyr::select(-cblAcc)
TmpTypeAcc <- "Cbl"
}
}
if (dim(ptsDirAcc)[1] > 1) {
PointTreeWIP <- ptsDirAcc %>% filter(type == "Tree")
if (verbose){
print(paste0("LoggedTrees : ",PointTreeWIP$origins[[1]]))
}
Lines[[l]] <- list("LineID" = NA,"LoggedTrees" = PointTreeWIP$origins, "TypeExpl" = TmpTypeAcc,"IdMachineZone" = ID_Access)
l <- l+1
}else{
#Compute Cost between points and Access points in cbl exploitation
CostDistEstCbl <- adjtopolcp(costSurface = CondSurf,
topography = DTMmean ,
pts = ptsWIPmax %>%
as_Spatial(),
slopeRdCond = SlopeCond,
paths = FALSE) [,1]
CostDistEstCbl <- CostDistEstCbl[(length(PointAcc$ID)+1):length(CostDistEstCbl)[1]]
CostDistEstCbl <- matrix(CostDistEstCbl,ncol = length(PointAcc$ID) )
#Attribute a least cost point access to each points
PointTree <- ptsWIPmax %>% filter(type != "Access" ) %>%
mutate(ptAccCbl = max.col(- CostDistEstCbl,ties.method = "first")) %>%
mutate(EstCostCbl = CostDistEstCbl)
if (winching == "2") {
CondSurfGrpl <- 1/CostRasterMeanGrpl
#Compute Cost between points and Access points in grpl exploitation
CostDistEstGrpl <- adjtopolcp(costSurface = CondSurfGrpl,
topography = DTMmean ,
pts = ptsWIPmax %>%
as_Spatial(),
slopeRdCond = SlopeCondGrpl,
paths = FALSE) [,1:length(PointAcc$ID)]
CostDistEstGrpl <- CostDistEstGrpl[(length(PointAcc$ID)+1):length(CostDistEstGrpl)[1]]
CostDistEstGrpl <- matrix(CostDistEstGrpl,ncol = length(PointAcc$ID) )
#Attribute a least cost point access to each points
PointTree <- PointTree %>%
mutate(ptAccGpl = max.col(-matrix(CostDistEstGrpl,ncol = length(PointAcc$ID) ),ties.method = "first")) %>%
mutate(EstCostGrpl = CostDistEstGrpl)
#Prioritize grpl exploitation if possible
CostDistEstGrpl[CostDistEstGrpl != Inf] <- 0
CostDistEstGrpl[CostDistEstGrpl == Inf] <- 1
for (j in 1:length(CostDistEstGrpl)) {
PointTree[j,"TypeAcc"] <- as.character(prod(CostDistEstGrpl[j,]))
}
PointTree$TypeAcc[PointTree$TypeAcc == "0"] <- "Grpl"
PointTree$TypeAcc[PointTree$TypeAcc == "1"] <- "Cbl"
PointTree$ptsAcc[PointTree$TypeAcc == "Grpl"] <- PointTree$ptAccGpl[PointTree$TypeAcc == "Grpl"]
PointTree$ptsAcc[PointTree$TypeAcc != "Grpl"] <- PointTree$ptAccCbl[PointTree$TypeAcc != "Grpl"]
PointTree$EstCost[PointTree$TypeAcc == "Grpl"] <- PointTree$EstCostGrpl[PointTree$TypeAcc == "Grpl"]
PointTree$EstCost[PointTree$TypeAcc != "Grpl"] <- PointTree$EstCostCbl[PointTree$TypeAcc != "Grpl"]
PointTree <- PointTree %>% arrange(desc(TypeAcc),EstCost)
#Define WIP point according to possible exploitation type
if (PointTree$TypeAcc[1] == "Grpl") {
#Reconstruct access points + selected tree in grpl exploitation
TmpPtsWIP <- ptsGrpl %>%
filter(IDpts == PointTree$IDpts[1]) %>%
st_union() %>%
st_cast("POINT", warn = FALSE) %>%
st_as_sf() %>%
mutate(type = "Overlay") %>%
mutate(ptsAcc = PointTree$ptsAcc[1]) %>%
mutate(IDpts = PointTree$IDpts[1]) %>%
mutate(origins = PointTree$origins[1])%>%
mutate(n.overlaps = PointTree$n.overlaps[1])
TmpTypeAcc <- "Grpl"
#Select adjacent grpl graph
SlopeCondRd <- SlopeCondGrpl
}else{
#if the selected point is not a grpl accessible point BUT not all grpl accessible point are done
#Shift to cbl exploitation
Grpl2CblFlag <- TRUE
#Reconstruct access points + selected tree in cbl exploitation
TmpPtsWIP <- ptsCbl %>%
filter(IDpts == PointTree$IDpts[1]) %>%
st_union() %>%
st_cast("POINT", warn = FALSE) %>%
st_as_sf() %>%
mutate(type = "Overlay") %>%
mutate(ptsAcc = PointTree$ptsAcc[1]) %>%
mutate(IDpts = PointTree$IDpts[1]) %>%
mutate(origins = PointTree$origins[1])%>%
mutate(n.overlaps = PointTree$n.overlaps[1])
#Select adjacent cbl graph
SlopeCondRd <- SlopeCond
TmpTypeAcc <- "Cbl"
}
}else{
#Reconstruct access points + selected tree in cbl exploitation
TmpPtsWIP <- ptsCbl %>%
filter(IDpts == PointTree$IDpts[1]) %>%
st_union() %>%
st_cast("POINT", warn = FALSE) %>%
st_as_sf() %>%
mutate(type = "Overlay") %>%
mutate(ptsAcc = PointTree$ptsAcc[1]) %>%
mutate(IDpts = PointTree$IDpts[1]) %>%
mutate(origins = PointTree$origins[1]) %>%
mutate(n.overlaps = PointTree$n.overlaps[1])
#Select adjacent cbl graph
SlopeCondRd <- SlopeCond
TmpTypeAcc <- "Cbl"
}
#filter WIP points in accessible scd trail area
# TmpPtsWIP <- st_set_crs(TmpPtsWIP,st_crs(ptsCbl)) # set crs from ptsCbl
TmpPtsWIP <- TmpPtsWIP %>%
filter(st_intersects(
st_geometry(TmpPtsWIP),st_geometry(AccessPolygons %>%
st_union()),sparse = FALSE))
#Reconstruct access points + selected points
TmpPtsWIP <- rbind( PointAcc %>%
st_union() %>%
st_cast("POINT", warn = FALSE)%>%
st_as_sf() %>%
mutate(type = "Access") %>%
mutate(ptsAcc = NA ) %>%
mutate(IDpts = NA) %>%
mutate(origins = NA)%>%
mutate(n.overlaps = NA),TmpPtsWIP)
#Compute Cost between all points and Access points
CostDistEstWIP <- adjtopolcp(costSurface = CondSurf,topography = DTMmean , pts = TmpPtsWIP %>% as_Spatial(),
slopeRdCond = SlopeCondRd,paths = FALSE)[,1]
CostDistEstWIP <- CostDistEstWIP[(length(PointAcc$ID)+1):length(CostDistEstWIP)[1]]
CostDistEstWIP <- matrix(CostDistEstWIP, ncol = length(PointAcc$ID))
PointTreeWIP <- TmpPtsWIP %>%
filter(type == "Overlay") %>%
mutate(EstCost = CostDistEstWIP) %>%
arrange(EstCost)
TmpPtsWIP <- rbind(TmpPtsWIP %>% filter(type == "Access") %>%
mutate(EstCost = NA),PointTreeWIP[1,])
TmpPathWIP <- adjtopolcp(costSurface = CondSurf,topography = DTMmean , pts = TmpPtsWIP %>% as_Spatial(),
slopeRdCond = SlopeCondRd,paths = FALSE)
TmpPathWIPCost <- TmpPathWIP[1:length(PointAcc$ID),length(PointAcc$ID)+1]
LCPathWIP <- max.col(t(-TmpPathWIPCost))
TmpPtsWIP <- rbind(TmpPtsWIP[LCPathWIP,],PointTreeWIP[1,])
TmpPathWIP <- adjtopolcp(costSurface = CondSurf,topography = DTMmean , pts = TmpPtsWIP %>% as_Spatial(),
slopeRdCond = SlopeCondRd,paths = TRUE)
if (TmpPathWIP[[1]][2,1] != 0) {
if (TmpTypeAcc == "Grpl") {
CostRasterMean <- rasterize(x = TmpPathWIP[[2]] ,
y = CostRasterMean ,
field = CostMatrix[[2]][[7]]$CostValue,update = TRUE)
CostRasterMeanGrpl <- rasterize(x = TmpPathWIP[[2]] ,
y = CostRasterMeanGrpl ,
field = CostMatrix[[2]][[7]]$CostValue,update = TRUE)
}else{
CostRasterMean <- rasterize(x = TmpPathWIP[[2]] ,
y = CostRasterMean ,
field = CostMatrix[[2]][[7]]$CostValue,update = TRUE)
}
}
pathLines[[k]] <- TmpPathWIP[[2]]
pathLines[[k]]@lines[[1]]@ID <- paste("Path",
"A",
LCPathWIP,
"NTree",
length(PointTreeWIP$origins[[1]]),
"T",
paste(as.character(unlist(PointTreeWIP$origins[[1]])),
collapse='-'),
sep = ".")
Lines[[l]] <- list("LineID" = k,"LoggedTrees" = PointTreeWIP$origins[[1]], "TypeExpl" = TmpTypeAcc,"IdMachineZone" = ID_Access)
k <- k +1
l <- l+1
pathLinesWIP[[ki]] <- TmpPathWIP[[2]]
pathLinesWIP[[ki]]@lines[[1]]@ID <- paste("Path",
"A",
LCPathWIP,
"NTree",
length(PointTreeWIP$origins[[1]]),
"T",
paste(as.character(unlist(PointTreeWIP$origins[[1]])),
collapse='-'),
sep = ".")
ki <- ki +1
pathsWIP <- do.call(rbind, pathLinesWIP)
pathsWIP <- pathsWIP %>% st_as_sf() %>% st_make_valid() %>%
filter(!st_is_empty(pathsWIP %>%
st_as_sf())) %>%
filter(st_geometry_type(.) %in% c("LINESTRING"))
if(dim(pathsWIP)[1]>0){
pathsWIP <- pathsWIP$geometry %>% as_Spatial()
if ( length(pathsWIP)==1) {
pathsWIP <- pathsWIP %>% st_as_sf() %>%
st_buffer(dist = advancedloggingparameters$ScndTrailWidth/2) %>% st_make_valid() %>%
st_union()
}else{
pathsWIP <- smoothtrails(paths = pathsWIP ,
plotmask = plotmask,
verbose = FALSE,
advancedloggingparameters = advancedloggingparameters)$SmoothedTrails %>%
st_union()
}
}else{
pathsWIP <- pathsWIP %>% st_union()
}
}
TreePts$Logged<- FALSE
for (j in 1:dim(TreePts)[1]) {
TreePts$Logged[j] <- any(PointTreeWIP$origins[[1]] %in% TreePts$idTree[j])
}
TreePts <- TreePts %>%
filter(Logged == FALSE) %>%
dplyr::select(-Logged)
if (dim(TreePts)[1]> 0) {
if (winching == "2") {
ptsGrpl <- TreePts %>% #def Grpl point
st_buffer(dist = advancedloggingparameters$GrappleLength) %>%
st_snap_to_grid(size = .2) %>% # avoid GEOS error (st-intersection issue)
#st_set_precision(1) %>%
st_temporaryintersection(topography = topography,
plotmask = plotmask,
advancedloggingparameters = advancedloggingparameters) %>%
st_make_valid() %>%
mutate(type = "Inter") %>%
dplyr::select(-ID)
for (OriginI in 1:length(ptsGrpl$origins)) {
for (OriginJ in 1:length(ptsGrpl$origins[[OriginI]])) {
IndexPts <- ptsGrpl$origins[[OriginI]][OriginJ]
ptsGrpl$origins[[OriginI]][OriginJ] <- TreePts$idTree[IndexPts]
}
}
ptsGrpl$isEmpty <- FALSE
for (h in 1:dim(ptsGrpl)[1]) {
if (ptsGrpl[h,] %>%
st_intersects(AccessMainTrails %>%
filter(ID == ID_Access),sparse = FALSE) == FALSE) {
ptsGrpl$isEmpty[h] <- TRUE
}else{
suppressWarnings(st_geometry(ptsGrpl[h,]) <- st_geometry(ptsGrpl[h,] %>%
st_intersection(AccessMainTrails %>%
filter(ID == ID_Access))))
}
}
#Filter polygons which intersect accessible area to second trails
ptsGrpl<- st_set_crs(ptsGrpl,st_crs(AccessPolygons)) # set crs from AccessPolygons
ptsGrpl <- ptsGrpl %>%
filter(isEmpty == FALSE) %>%
dplyr::select(-isEmpty) %>%
mutate(IDpts = paste0("I.", row_number()))
ptsWIP <- ptsGrpl %>% #def Grpl point as WIP point
st_set_agr("constant") %>%
st_centroid()
ptsCbl <- TreePts %>% #def cbl polygons
st_buffer(dist = advancedloggingparameters$CableLength) %>%
st_snap_to_grid(size = 1) %>%# avoid GEOS error (st-intersection issue)
st_set_precision(1) %>%
st_temporaryintersection(topography = topography,
plotmask = plotmask,
advancedloggingparameters = advancedloggingparameters) %>%
st_make_valid() %>%
mutate(type = "Inter")%>%
dplyr::select(-ID)
for (OriginI in 1:length(ptsCbl$origins)) {
for (OriginJ in 1:length(ptsCbl$origins[[OriginI]])) {
IndexPts <- ptsCbl$origins[[OriginI]][OriginJ]
ptsCbl$origins[[OriginI]][OriginJ] <- TreePts$idTree[IndexPts]
}
}
ptsCbl$isEmpty <- FALSE
for (h in 1:dim(ptsCbl)[1]) {
if (ptsCbl[h,] %>%
st_intersects(st_geometry(AccessMainTrails %>%
filter(ID == ID_Access) %>%
st_union()),sparse = FALSE) == FALSE) {
ptsCbl$isEmpty[h] <- TRUE
}else{
suppressWarnings(st_geometry(ptsCbl[h,]) <- st_geometry(ptsCbl[h,] %>%
st_intersection(st_geometry(AccessMainTrails %>%
filter(ID == ID_Access)))))
}
}
ptsCbl <- st_set_crs(ptsCbl,st_crs(AccessPolygons)) #set crs from AccessPolygons
ptsCbl <- ptsCbl %>%
filter(isEmpty == FALSE) %>%
dplyr::select(-isEmpty) %>%
mutate(IDpts = paste0("I.",row_number())) #Filter polygons which intersect accessible area to second trails
ptsWIPCbl <- ptsCbl %>%#Convert polygons to centroid
st_set_agr("constant") %>%
st_centroid()
ptsWIP <- ptsWIP %>%
arrange(desc(n.overlaps))
#Select adjacent grpl graph
SlopeCondRd <- SlopeCondGrpl
}else{
ptsCbl <- TreePts %>% #def cbl point
st_buffer(dist = advancedloggingparameters$CableLength) %>%
st_snap_to_grid(size = 1) %>% # avoid GEOS error (st-intersection issue)
st_set_precision(1) %>%
st_temporaryintersection(topography = topography,
plotmask = plotmask,
advancedloggingparameters = advancedloggingparameters) %>%
st_make_valid() %>%
mutate(IDpts = paste0("I.",row_number())) %>%
mutate(type = "Inter") %>% dplyr::select(-ID)
ptsCbl <- st_set_crs(ptsCbl,st_crs(AccessPolygons)) #set crs from AccessPolygons
for (OriginI in 1:length(ptsCbl$origins)) {
for (OriginJ in 1:length(ptsCbl$origins[[OriginI]])) {
IndexPts <- ptsCbl$origins[[OriginI]][OriginJ]
ptsCbl$origins[[OriginI]][OriginJ] <- TreePts$idTree[IndexPts]
}
}
ptsCbl$isEmpty <- FALSE
for (h in 1:dim(ptsCbl)[1]) {
if (ptsCbl[h,] %>%
st_intersects(AccessMainTrails %>%
filter(ID == ID_Access),sparse = FALSE) == FALSE) {
ptsCbl$isEmpty[h] <- TRUE
}else{
suppressWarnings(st_geometry(ptsCbl[h,]) <- st_geometry(ptsCbl[h,] %>%
st_intersection(AccessMainTrails %>%
filter(ID == ID_Access))))
}
}
ptsWIP <- ptsCbl %>%
filter(isEmpty == FALSE) %>%
dplyr::select(-isEmpty) %>%#def cbl point as WIP point
st_set_agr("constant") %>%
st_centroid()
ptsWIP <- ptsWIP %>% #filter cbl intersection centroid point out plot
filter(st_intersects(st_geometry(ptsWIP),
st_geometry(plotmask %>% st_as_sf()),
sparse = FALSE)) %>%
arrange(desc(n.overlaps))
#Select adjacent cbl graph
SlopeCondRd <- SlopeCond
}
}
RemainTree <- dim(TreePts)[1]
}
if (RemainTree > 0L){
stop("Failed to compute 2nd trails")
}
}
}
MainTrailsAccess <- AccessPointAll %>%
st_set_crs(st_crs(topography))
inventory <- treeselectionoutputs$inventory
# initial inventory
inventory0 <- inventory
lines <- do.call(rbind, Lines)
# WinchingMachine info in the inventory
TrailsIdentity_df <- as.data.frame(lines) %>%
unnest(cols = c(LineID, LoggedTrees, TypeExpl, IdMachineZone)) # unnesting flattens it back out into regular columns
TrailsIdentity_df <- TrailsIdentity_df %>%
dplyr::select(LoggedTrees, TypeExpl, IdMachineZone) %>%
rename(idTree = LoggedTrees) %>%
rename(WinchingMachine = TypeExpl) %>%
unique()
if(any(duplicated(TrailsIdentity_df$idTree))) # trees logged several times?
stop("Error: duplicated tree in logged trees")
inventory <- inventory %>%
left_join(TrailsIdentity_df, by = "idTree")
paths <- do.call(rbind, pathLines)
paths <- paths %>% st_as_sf() %>% st_make_valid() %>%
filter(!st_is_empty(paths %>% st_as_sf())) %>%
filter(st_geometry_type(.) %in% c("LINESTRING"))
if (dim(paths)[1]>0) {
paths <- paths$geometry %>% as_Spatial()
if ( length(paths)==1) {
SmoothedTrails <- paths %>% st_as_sf() %>%
st_buffer(dist = advancedloggingparameters$ScndTrailWidth/2) %>% st_make_valid() %>%
st_union()%>% st_set_crs(st_crs(topography))
TrailsDensity <- (SmoothedTrails %>%
st_intersection(plotmask %>% st_as_sf()) %>%
st_area / advancedloggingparameters$ScndTrailWidth)/(plotmask %>% st_as_sf() %>% st_area() /10000)
}else{
secondtrails <- smoothtrails(paths,
plotmask,
verbose = verbose,
advancedloggingparameters = advancedloggingparameters)
SmoothedTrails <- secondtrails$SmoothedTrails %>% st_set_crs(st_crs(topography))
TrailsDensity <- secondtrails$TrailsDensity
}
#
#
# # Records the dead trees
if(fuel == "0"){
if (!("DeathCause" %in% names(inventory))){
inventory <- inventory %>%
add_column(DeathCause = NA) # if "DeathCause" column doesnt exist create it
}
DeadTrees <- suppressWarnings(st_intersection(
st_make_valid(st_set_crs(st_as_sf(inventory, coords = c("Xutm", "Yutm")), st_crs(topography))),
st_make_valid(st_as_sf(SmoothedTrails) %>% st_union()) # "make valid" to avoid self-intersection
)) %>%
add_column(DeadTrees = "1") %>%
dplyr::select(idTree, DeadTrees) %>% st_drop_geometry() # remove geometry column (sf to data.frame)
DeadTrees <- unique(DeadTrees)
inventory <- inventory %>%
left_join(DeadTrees, by = "idTree") %>%
mutate(DeathCause = ifelse(is.na(DeathCause) & Selected != "1" & DeadTrees == "1",
"2ndtrail", DeathCause)) %>% # Damage trees
dplyr::select(-DeadTrees)
}
}else{
paths <- NULL
SmoothedTrails <- MainTrailsAccess %>% st_union() %>%
st_buffer(dist = advancedloggingparameters$MaxMainTrailWidth/2) %>% st_set_crs(st_crs(topography))
TrailsDensity <- 0
if(fuel == "0"){
if (!("DeathCause" %in% names(inventory))){
inventory <- inventory %>%
add_column(DeathCause = NA) # if "DeathCause" column doesnt exist create it
}
DeadTrees <- suppressWarnings(sf::st_intersection(
sf::st_make_valid(st_set_crs(st_as_sf(inventory, coords = c("Xutm", "Yutm")), st_crs(topography))),
sf::st_make_valid(st_as_sf(SmoothedTrails) %>% st_union()) # "make valid" to avoid self-intersection
)) %>%
add_column(DeadTrees = "1") %>%
dplyr::select(idTree, DeadTrees) %>% st_drop_geometry()
DeadTrees <- unique(DeadTrees)
inventory <- inventory %>%
left_join(DeadTrees, by = "idTree") %>%
mutate(DeathCause = ifelse(is.na(DeathCause) & Selected != "1" & DeadTrees == "1",
"2ndtrail", DeathCause)) %>% # Damage trees
dplyr::select(-DeadTrees)
}
}
if(nrow(inventory) != nrow(inventory0))
stop("The number of rows between the input inventory and the output inventory
of the function secondtrailsopening() is not the same.The function must be corrected.")
# OUTPUTS
if (WinchingInit == "2") {
secondtrails <- list("inventory" = inventory,
"SmoothedTrails" = SmoothedTrails,
"TrailsDensity" = TrailsDensity,
"TrailsIdentity" = lines, # "LineID","LoggedTrees", "TypeExpl"
"MainTrailsAccess" = MainTrailsAccess,
"RawSecondTrails" = paths,
"CostRasterAgg" = list("CostRasterMean" = CostRasterMean,"CostRasterMeanGrpl" = CostRasterMeanGrpl))
}else{
secondtrails <- list("inventory" = inventory,
"SmoothedTrails" = SmoothedTrails,
"TrailsDensity" = TrailsDensity,
"TrailsIdentity" = lines, # "LineID","LoggedTrees", "TypeExpl"
"MainTrailsAccess" = MainTrailsAccess,
"RawSecondTrails" = paths,
"CostRasterAgg" = list("CostRasterMean" = CostRasterMean,"CostRasterMeanGrpl" = NULL))
}
return(secondtrails)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.