R/snap_tools.R

In Jean-Romain/ALSroads: Road alignment and measure from airborne LiDAR data for MFFP

#' Snap road endings
#'
#' Snap road endings together, with or without prior knowledge on road connections, to post-process
#' and correct inaccuracies generated by \link{measure_roads} and ensure getting topologically valid
#' network. If argument 'ref' is missing the method is very basic. The method using a reference map
#' is more advanced.
#'
#' @param roads  multiple lines (\code{sf} format). Corrected but unconnected roads.
#' @param ref  multiple lines (\code{sf} format). Original non-corrected but connected roads. Can be
#' missing.
#' @param tolerance numeric (distance unit). Tolerance value used to snap the road endings.
#' @param field  character. Unique identifier field in both road datasets. Relevant only if 'ref' is
#' not missing.
#' @param updatable logical. Vector identifing which road segments need to be snapped (default is all
#' segments). Relevant only if 'ref' is not missing.
#'
#' @return Named list with \code{roads} being the same object as \code{roads} but with corrected endings
#' such that roads are connected and \code{warnings} being either \code{NULL} or a \code{sf POINT}
#' identifying nodes with issues.
#' @examples
#' f <- system.file("extdata", "j53e_network.gpkg", package="ALSroads")
#'
#' ref <- sf::st_read(f, layer = "original")  # input of measure_roads
#' cor <- sf::st_read(f, layer = "corrected") # output of measure_roads
#' res <- st_snap_lines(cor, ref, field = "OBJECTID")
#'
#' plot(sf::st_geometry(ref), xlim = c(260800, 261000), ylim = c(5250400, 5250650), col = "red")
#' plot(sf::st_geometry(cor), col = "blue", add = TRUE)
#' plot(sf::st_geometry(res$roads), col = "green", add = TRUE, lwd = 2)
#'
#' domain <- "https://servicesmatriciels.mern.gouv.qc.ca:443"
#' path <- "/erdas-iws/ogc/wmts/Inventaire_Ecoforestier/Inventaire_Ecoforestier/default/"
#' tiles <- "GoogleMapsCompatibleExt2:epsg:3857/{z}/{y}/{x}.jpg"
#' url <- paste0(domain, path, tiles)
#' m = mapview::mapview(list(ref, cor, res$roads),
#'   layer.name = c("Inaccurate", "Corrected", "Snapped"),
#'   color = c("red", "blue", "green"), map.type = "Esri.WorldImagery")
#' leaflet::addTiles(m@map, url)
#' @export
st_snap_lines = function(roads, ref, tolerance = 30, field = NULL, updatable = NULL)
{
  if (missing(ref))
    return(simple_snap(roads, tolerance))

  if (!missing(ref) & is.null(field))
    stop("Argument 'field' cannot be NULL if a reference topology is provided", call. = FALSE)

  if (is.null(updatable))
    updatable <- rep(TRUE, nrow(roads))

  return(advanced_snap(roads, ref, field, tolerance, updatable))
}

simple_snap <- function(roads, tolerance)
{
  pts_warn <- NULL
  end   <- lwgeom::st_endpoint(roads)
  start <- lwgeom::st_startpoint(roads)
  ends  <- c(start, end)

  u <- sf::st_is_within_distance(ends, ends, tolerance)
  u <- lapply(u, sort)
  u <- Filter(function(x) { length(x) > 1 }, u)
  u <- unique(u)
  u <- lapply(u, function(x) { ends[x] })
  u <- lapply(u, function(x) { sf::st_sfc(sf::st_point(colMeans(sf::st_coordinates(x)))) })
  u <- do.call(c, u)

  if (is.null(u))
  {
    dist_unit <- sf::st_crs(roads)$units
    warning(glue::glue("No roads could be snapped with the tolerance used ({tolerance} {dist_unit}). Original roads returned."), call. = FALSE)
    return(list(roads = roads, warnings = NULL))
  }

  sf::st_crs(u) <- sf::st_crs(roads)

  v <- sf::st_is_within_distance(u, start, tolerance)
  for (i in seq_along(v))
  {
    ids <- v[[i]]
    for(j in ids)
    {
      sf::st_geometry(roads[j,])[[1]][1, 1] <- u[i][[1]][[1]]
      sf::st_geometry(roads[j,])[[1]][1, 2] <- u[i][[1]][[2]]
    }
  }

  w <- sf::st_is_within_distance(u, end, tolerance)
  for (i in seq_along(w))
  {
    ids <- w[[i]]
    for(j in ids)
    {
      n <- nrow(sf::st_geometry(roads[j,])[[1]])
      sf::st_geometry(roads[j,])[[1]][n, 1] <- u[i][[1]][[1]]
      sf::st_geometry(roads[j,])[[1]][n, 2] <- u[i][[1]][[2]]
    }
  }

  end   <- lwgeom::st_endpoint(roads)
  start <- lwgeom::st_startpoint(roads)
  idx_idem <- which(end == start)

  if (length(idx_idem))
  {
    warning("Some roads had both of their ends being snapped together.", call. = FALSE)
    pts_warn <- sf::st_sf(start[idx_idem], sf::st_drop_geometry(roads[idx_idem,]))
  }

  return(list(roads = roads, warnings = pts_warn))
}

advanced_snap <- function(roads, ref, field, tolerance, updatable)
{
  X <- Y <- id <- SCORE <- CLASS <- UPDATABLE <- NULL

  IDs1 <- sort(unique(roads[[field]]))
  IDs2 <- sort(unique(ref[[field]]))
  if (length(IDs1) != length(roads[[field]])) stop("Values in unique identifier field are not unique for 'roads'.", call. = FALSE)
  if (length(IDs2) != length(ref[[field]])) stop("Values in unique identifier field are not unique for 'ref'.", call. = FALSE)
  if (!all(IDs1 == IDs2)) stop("Values in unique identifier field are not the same in both road datasets.", call. = FALSE)
  if (!"SCORE" %in% colnames(roads)) stop("'SCORE' column is missing from 'roads'.", call. = FALSE)
  if (!"CLASS" %in% colnames(roads)) stop("'CLASS' column is missing from 'roads'.", call. = FALSE)
  if (!all(is.logical(updatable)) | length(updatable) != nrow(roads)) stop("'updatable' must be a logical vector of the same length of 'roads'.", call. = FALSE)

  dist_unit <- sf::st_crs(roads)$units
  df_pts_warn <- data.frame(X = as.numeric(), Y = as.numeric(), message = as.character(), IDs = as.character())
  pts_warn <- NULL

  # Arrange both road datasets to make sure
  # that line indices will match between them
  roads <- roads |>
    dplyr::mutate(UPDATABLE = updatable) |>
    dplyr::arrange(field)
  ref <- dplyr::arrange(ref, field)

  prepare_data <- function(roads, end = FALSE)
  {
    if (end) {
      get_end <- lwgeom::st_endpoint
    } else {
      get_end <- lwgeom::st_startpoint
    }

    out <- get_end(roads) |>
      sf::st_coordinates() |>
      dplyr::as_tibble() |>
      dplyr::mutate(pos = 1:nrow(roads)) |>
      dplyr::mutate(end = end) |>
      dplyr::mutate(id = paste0(roads[[field]], "_", end))

    return(out)
  }

  # Extraction of roads end points coordinates along with heading and
  # association with a unique identifier

  # Corrected roads
  tb_endpoint <- prepare_data(roads, TRUE)
  tb_startpoint <- prepare_data(roads, FALSE)

  tb_ends_roads <- rbind(tb_endpoint, tb_startpoint) |>
    dplyr::mutate(SCORE = rep(roads[["SCORE"]], 2)) |>
    dplyr::mutate(CLASS = rep(roads[["CLASS"]], 2)) |>
    dplyr::mutate(UPDATABLE = rep(roads[["UPDATABLE"]], 2))

  # Non-corrected roads
  tb_endpoint   <- prepare_data(ref, TRUE)
  tb_startpoint <- prepare_data(ref, FALSE)

  ls_heading <- lapply(sf::st_geometry(ref), st_ends_heading)
  heading_tail_head <- c(sapply(ls_heading, utils::tail, 1),
                         sapply(ls_heading, utils::head, 1))

  tb_ends_ref <- rbind(tb_endpoint, tb_startpoint) |>
    dplyr::mutate(heading = heading_tail_head)


  # List of distinct connected nodes in the non-corrected/reference road dataset
  tb_nodes_grouped <- dplyr::group_by(tb_ends_ref, X, Y)

  if (nrow(tb_nodes_grouped) == 0)
  {
    warning("No road junction has been found in reference road dataset. Original roads returned.", call. = FALSE)
    return(list(roads = roads, warnings = NULL))
  }

  ls_tb_nodes <- tb_nodes_grouped |>
    dplyr::group_split() |>
    lapply(function(x) if (nrow(x) > 1) x)
  ls_tb_nodes <- ls_tb_nodes[!sapply(ls_tb_nodes, is.null)]


  # Try snapping together each segment at each node found.
  # This loop can't be vectorised because of the geometry
  # change occuring at both ends of a same road segment
  for (tb_node in ls_tb_nodes)
  {
    IDs <- tb_node[["id"]]
    tb_node_ref <- dplyr::filter(tb_ends_ref, id %in% IDs)
    tb_node_cor <- dplyr::filter(tb_ends_roads, id %in% IDs)

    # Prepare IDs string for potential warning message
    pos <- unique(tb_node[["pos"]])
    IDs_node <- roads[pos,][[field]]
    IDs_glued <- glue::glue_collapse(IDs_node, ", ")


    # Will only try snapping if at least one of the segments
    # in the group has been corrected as they will already
    # be connected if none has been corrected
    been_corrected <- any(tb_node_cor[["UPDATABLE"]])
    if (!been_corrected) next


    # Find the two segments that are the most
    # likely to be the extension of each other
    # Those two are marked out as the "bridge"
    bridge_ids <- tb_node_ref |>
      dplyr::left_join(dplyr::select(tb_node_cor, id, SCORE, CLASS), by = "id") |>
      find_best_connexion()

    tb_node_bridge <- dplyr::filter(tb_node_cor, id %in% bridge_ids)

    # Make sure the gap between the two segments that
    # will make the bridge isn't exceeding tolerance otherwise
    # skip snapping this node and throw warning
    gap <- sqrt(diff(tb_node_bridge[["X"]])^2 + diff(tb_node_bridge[["Y"]])^2)
    if (gap/2 > tolerance)
    {
      warning(glue::glue("Impossible to connect together roads with '{field}' {IDs_glued}; distance from expected junction exceeds tolerance ({round(gap/2,1)} > {tolerance} {dist_unit})."), call.=FALSE)
      df_pts_warn <- data.frame(tb_node_ref[1, c("X","Y")], message = "Distance from expected junction exceeds tolerance", IDs = IDs_glued) |>
        rbind(df_pts_warn)
      next
    }


    # If there is only two segments, there is no need
    # to try finding a better junction point with a
    # non-existent third or fourth segment...
    if (length(IDs) == 2)
    {
      # Update geometry of the two corrected segments composing the bridge
      for (j in 1:2)
      {
        pos <- tb_node_bridge[j,][["pos"]]
        n <- ifelse(tb_node_bridge[j,][["end"]], nrow(sf::st_coordinates(roads[pos,])), 1)
        sf::st_geometry(roads[pos,])[[1]][n, 1] <- mean(tb_node_bridge[["X"]])
        sf::st_geometry(roads[pos,])[[1]][n, 2] <- mean(tb_node_bridge[["Y"]])
      }
      next
    }


    # If the node contains a loop, special handling is needed
    tb_pos <- table(tb_node[["pos"]])
    if (any(tb_pos > 1))
    {
      if (length(tb_pos) == 2)
      {
        # Update geometry in the case of one loop and one normal road
        for (j in 1:3)
        {
          pos <- tb_node[j,][["pos"]]
          n <- ifelse(tb_node[j,][["end"]], nrow(sf::st_coordinates(roads[pos,])), 1)
          sf::st_geometry(roads[pos,])[[1]][n, 1] <- mean(tb_node_bridge[["X"]])
          sf::st_geometry(roads[pos,])[[1]][n, 2] <- mean(tb_node_bridge[["Y"]])
        }
      }
      else
      {
        # Raise warning and skip in the case of one or more loops and one or more normal road
        # Could be adressed but this case will be so rare that I don't think
        # it justifies the complexity that would need to be added to the rest of the script
        # One way of doing it would be to split each loop before searching for the best
        # junction point and then only at the end recombining each loop.
        warning(glue::glue("Impossible to connect together roads with '{field}' {IDs_glued}; junction contains one loop and at least two other roads."), call.=FALSE)
        df_pts_warn <- data.frame(tb_node_ref[1, c("X","Y")], message = "Junction contains one loop and at least two other roads", IDs = IDs_glued) |>
          rbind(df_pts_warn)
      }
      next
    }


    # Merge the two segments forming the bridge.
    # This will allow for an easier split when the
    # exact position of the junction will be found.
    bridge <- make_bridge(roads, tb_node_bridge)


    # Find intersection points between the remaining
    # segments and the bridge. The intersection point
    # is given as the distance on the path between
    # the beginning of the bridge and the intersection point.
    remain_ids <- IDs[!(IDs %in% bridge_ids)]
    tb_node_remain <- dplyr::filter(tb_node_cor, id %in% remain_ids)

    pos <- tb_node_remain[["pos"]]
    end <- tb_node_remain[["end"]]
    head_tail <- sapply(as.numeric(end) + 1, function(x) c("HEAD", "TAIL")[x])

    lines_extended <- roads[pos,] |>
        sf::st_geometry() |>
        mapply(999999, head_tail, FUN = st_extend_line, SIMPLIFY = FALSE) |>
        sf::st_sfc(crs = sf::st_crs(roads))

    dist_bridge <- sapply(lines_extended,
                          distance_line_intersection,
                          bridge[["bridge"]],
                          bridge[["junction"]])


    # Make sure all intersections occur within the tolerance
    # value of the mean junction otherwise skip snapping
    # this node and throw warning
    dist_junction_bridge <- st_split_at_point(bridge[["junction"]], bridge[["bridge"]])[1] |> sf::st_length() |> as.numeric()

    no_intersection <- sapply(dist_bridge, is.na)
    dist_junction_diff <- abs(dist_bridge - dist_junction_bridge)
    dist_junction_diff[no_intersection] <- Inf
    dist_max <- max(dist_junction_diff)

    if (dist_max > tolerance)
    {
      warning(glue::glue("Impossible to connect together roads with '{field}' {IDs_glued}; distance from expected junction exceeds tolerance ({round(dist_max,1)} > {tolerance} {dist_unit})."), call.=FALSE)
      df_pts_warn <- data.frame(tb_node_ref[1, c("X","Y")], message = "Distance from expected junction exceeds tolerance", IDs = IDs_glued) |>
        rbind(df_pts_warn)
      next
    }

    # Find coordinates of the point where
    # all remaining segments intersect the bridge
    # based on the mean intersection distance
    # along the bridge from its begining
    dist_junction_mean <- mean(dist_bridge)
    mean_junction <- st_point_on_line(dist_junction_mean, bridge[["bridge"]])

    # Make sure that the mean intersection point doesn't
    # occur at the end of the bridge. Really rare, but
    # can happen in case of a short by-pass road misplaced
    # but already connected at one end of the bridge
    if (dist_junction_mean == 0)
    {
      warning(glue::glue("Impossible to connect together roads with '{field}' {IDs_glued}; junction weirdly occurs at exactly one end of a road."), call.=FALSE)
      df_pts_warn <- data.frame(tb_node_ref[1, c("X","Y")], message = "Junction weirdly occurs at exactly one end of a road", IDs = IDs_glued) |>
        rbind(df_pts_warn)
      next
    }

    # Split bridge geometry to recreate the original
    # two segments composing it and get coordinates
    # of the final junction point
    bridge_geometries <- st_split_at_point(mean_junction, bridge[["bridge"]])

    if (length(bridge_geometries) != 2)
    {
      warning(glue::glue("Impossible to connect together roads with '{field}' {IDs_glued}; issue occured during a splitting operation."), call.=FALSE)
      df_pts_warn <- data.frame(tb_node_ref[1, c("X","Y")], message = "Issue occured during a splitting operation", IDs = IDs_glued) |>
        rbind(df_pts_warn)
      next
    }

    coords_junction <- sf::st_coordinates(bridge_geometries[2])[1,-3]

    # Reorder vertices in order to match the original
    # two segments forming the bridge
    for (k in which(bridge[["reversed"]]))
    {
      bridge_geometries[k][[1]] <- sf::st_coordinates(bridge_geometries[k])[,-3] |>
        apply(2, rev) |>
        sf::st_linestring()
    }

    # Update geometry of the two corrected segments composing the bridge
    pos <- tb_node_bridge[["pos"]]
    sf::st_geometry(roads[pos,]) <- bridge_geometries


    # Update coordinates of remaining segments
    # by splitting in order to avoid potential
    # overlapping segment if it was already crossing
    # the bridge before elongation
    idx_valid_intersection <- which(!no_intersection)
    for (j in idx_valid_intersection)
    {
      pos <- tb_node_remain[j,][["pos"]]
      end <- tb_node_remain[j,][["end"]]
      head_tail <- c("HEAD", "TAIL")[as.numeric(end) + 1]

      road_extended <- roads[pos,] |>
        sf::st_geometry() |>
        st_extend_line(999999, head_tail)

      road_split <- road_extended |>
        lwgeom::st_split(bridge[["bridge"]]) |>
        sf::st_collection_extract("LINESTRING")


      # Find vertices coordinates corresponding
      # to the bulk of the original segment which
      # aren't those between the extended vertex
      # and the junction vertex with the bridge
      coords_split <- sf::st_coordinates(road_split)
      coords_extended <- sf::st_coordinates(road_extended)
      n <- ifelse(end, nrow(coords_extended), 1)
      x_match_extended <- coords_split[,"X"] == coords_extended[n, "X"]
      y_match_extended <- coords_split[,"Y"] == coords_extended[n, "Y"]
      idx_match_extended <- which(x_match_extended & y_match_extended)

      # Check if a split really occured as it might
      # not in the weird and rare case where the segment is near
      # enough of the mean jonction but at the same time
      # doesn't cross the bridge like if one end of the segment
      # is exactly at the junction point (I saw this only once...)
      if (length(road_split) == 1)
      {
        ID_problem <- roads[pos,][[field]]
        warning(glue::glue("Impossible to connect road with '{field}' {ID_problem} inside node of {IDs_glued}; the road doesn't intersect other roads but weirdly is still very close."), call.=FALSE)
        df_pts_warn <- data.frame(tb_node_ref[1, c("X","Y")], message = "Road doesn't seem to intersect with others inside node", IDs = ID_problem) |>
          rbind(df_pts_warn)
        next
      }


      # Some tolerance must be added due to slight imprecision
      # in the coordinates returned by st_point_on_line()
      coords_intersect <- dist_bridge[j] |>
        st_point_on_line(bridge[["bridge"]]) |>
        sf::st_coordinates()

      x_match_junc <- (coords_split[,"X"] > (coords_intersect[1,"X"] - 0.01)) &
                      (coords_split[,"X"] < (coords_intersect[1,"X"] + 0.01))
      y_match_junc <- (coords_split[,"Y"] > (coords_intersect[1,"Y"] - 0.01)) &
                      (coords_split[,"Y"] < (coords_intersect[1,"Y"] + 0.01))
      idx_match_junc <- which(x_match_junc & y_match_junc)

      # Remove all coordinates between the extended vertex (included)
      # and the junction vertex (excluded)
      idx_near_junc <- ifelse(idx_match_junc[1] > idx_match_extended, idx_match_junc[1], idx_match_junc[2])
      coords_keep <- coords_split[-(idx_near_junc:idx_match_extended),-3]

      # Edit the coordinates at (mean) splitting point
      idx_junction <- ifelse(end, nrow(coords_keep), 1)
      coords_keep[idx_junction,] <- coords_junction

      # Remove duplicates vertices created by the splitting operation
      idx_duplicated <- which((diff(coords_keep[,"X"]) == 0) & (diff(coords_keep[,"Y"]) == 0))
      if (length(idx_duplicated)) coords_keep <- coords_keep[-c(idx_duplicated, idx_duplicated+1),]


      # Update geometry of the corrected segment
      sf::st_geometry(roads[pos,]) <- coords_keep |>
        sf::st_linestring() |>
        sf::st_sfc(crs = sf::st_crs(roads))
    }
  }

  if (nrow(df_pts_warn))
    pts_warn <- sf::st_as_sf(x = df_pts_warn, coords = c("X","Y"), crs = sf::st_crs(roads))

  return(list(roads = dplyr::select(roads, -UPDATABLE), warnings = pts_warn))
}


#' Find intersection point between two lines as distance
#'
#' Find the closest intersection point between two lines
#' relative to a reference point on the second line.
#'
#' @param crossing_line  line (\code{sfc} format)
#' @param reference_line  line (\code{sfc} format). Reference line against which \code{crossing_line} intersection will be checked.
#' @param reference_point  point (\code{sfc} format). Point that must be on \code{reference_line} and that will be used to find
#' the closest intersection point of \code{crossing_line} to itself.
#'
#' @return Distance (\code{numeric}) between the beginning of \code{reference_line} and the closest intersection point to \code{reference_point}
#' on its path. \code{NA} is returned if no intersection occur.
#' @noRd
distance_line_intersection <- function(crossing_line, reference_line, reference_point)
{
  # Removing CRS information of reference_line/point allow
  # using distance_line_intersection() with functions
  # of the apply() family. Those function subset the
  # object passed as first argument meaning that if
  # a sfc_LINESTRING object is provided, only a
  # LINESTRING object will get passed to the function
  # which will cause an error later at sf::st_intersection()
  sf::st_crs(reference_line) <- NA
  sf::st_crs(reference_point) <- NA

  split_at_crossing <- lwgeom::st_split(reference_line, crossing_line) |> sf::st_collection_extract("LINESTRING")
  
  if (length(split_at_crossing) == 1) return(NA)
  
  # Compute distance between each crossing of the
  # reference line and the reference point on this line
  # and only keep the closest crossing
  dist_intersect <- sapply(head(split_at_crossing, -1), sf::st_length) |> cumsum()
  dist_point <- st_split_at_point(reference_point, reference_line)[1] |> sf::st_length()
  dist_intersect_closest <- dist_intersect[which.min(abs(dist_intersect - dist_point))]

  return(dist_intersect_closest)
}


#' Find best connexion between two lines at a junction
#'
#' At a junction with multiple lines, find which two lines are the most likely
#' extension of each other. The prime factor is the angle formed between
#' each pair of lines (the flatter the better).
#'
#' @param tb_node  \code{tibble} containing info about the heading, score and class of all lines
#' at a same junction.
#'
#' @return IDs (\code{character}) of the pair of lines having the best fit.
#' @noRd
find_best_connexion <- function(tb_node)
{
  CLASS <- NULL

  # Update CLASS 0 to CLASS 5 in order to get a constant
  # quality gradient from 1 to 5 for an easier sort
  tb_node <- dplyr::mutate(tb_node, CLASS = ifelse(CLASS == 0, 5, CLASS))

  # Produce lines permutations (as row number)
  m_row <- utils::combn(1:nrow(tb_node), 2)

  # Compute minimal angle difference between each lines
  min_diff_angle <- function(x, y)
  {
    a <- x - y
    if (a > 180) a <- 360 - a
    if (a < 0) a <- 360 + a
    min(a, 360 - a)
  }
  comb_angle <- apply(m_row, 2, function(x) min_diff_angle(tb_node[["heading"]][x[1]], tb_node[["heading"]][x[2]]))

  # Give a chance to lines that have an angle close to "flattest" one to
  # be considered "as flat". The flatness of the angle being the the prime
  # driver of the selection, we want to avoid a CLASS 0 line taking
  # precedence over a CLASS 1 line for a mere 5° difference
  tolerance_angle <- 7.5  # Allow angle up to 7.5° more
  idx_max <- which.max(comb_angle)[1]
  comb_angle[comb_angle >= comb_angle[idx_max] - tolerance_angle] <- 180

  comb_class <- comb_score <- NULL

  # Build similarity table using all permutations
  tb_similarity <-
    dplyr::tibble(
      pairs        = 1:ncol(m_row),
      comb_angle   = comb_angle,
      comb_class   = apply(m_row, 2, function(x) sum(tb_node[["CLASS"]][x]) ),
      comb_score   = apply(m_row, 2, function(x) sum(tb_node[["SCORE"]][x], na.rm = TRUE) )) |>
    dplyr::arrange(dplyr::desc(comb_angle), comb_class, dplyr::desc(comb_score))

  # The best fit is based on (in order):
  #  - the angle formed between the segments (the flatter the better)
  #  - the class (corrected roads are favored)
  #  - the score (best overal score is favored)
  idx_best <- tb_similarity[1,][["pairs"]]
  ids_best <- tb_node[m_row[,idx_best],][["id"]]

  return(ids_best)
}


#' Combine two connected lines
#'
#' Find the mean coordinates of the two closest ends and combine
#' the two lines at this point to create a single continuous line
#' while maintaining history of the vertex order in the two
#' original lines.
#'
#' @param roads  lines (\code{sf} format). Corrected roads.
#' @param tb_node_bridge  \code{tibble} containing info about the two lines to be connected. Must contain
#' the corresponding row index in \code{roads} (field 'pos') as well as boolean value indicating if the vertex to connect correspond
#' to the end/tail one (field 'end').
#'
#' @return Named list containing the created \code{bridge} (as \code{sfc_LINESTRING}), the \code{junction} point (as \code{sfc_POINT})
#' and a vector indicating which line has been \code{reversed}.
#' @noRd
make_bridge <- function(roads, tb_node_bridge)
{
  end_1 <- tb_node_bridge[1,][["end"]]
  end_2 <- tb_node_bridge[2,][["end"]]
  pos_1 <- tb_node_bridge[1,][["pos"]]
  pos_2 <- tb_node_bridge[2,][["pos"]]
  coord_1 <- sf::st_coordinates(roads[pos_1,])[,-3]
  coord_2 <- sf::st_coordinates(roads[pos_2,])[,-3]
  rev_1 <- FALSE
  rev_2 <- FALSE

  # Sometimes, reversing direction of one of the
  # two segments is needed in order to obtain a
  # "continuous flow" of the vertices
  if (!end_1)
  {
    rev_1 <- TRUE
    coord_1 <- apply(coord_1, 2, rev)
    if (end_2)
    {
      rev_2 <- TRUE
      coord_2 <- apply(coord_2, 2, rev)
    }
  } else {
    if (end_2)
    {
      rev_2 <- TRUE
      coord_2 <- apply(coord_2, 2, rev)
    }
  }

  # Mean coordinates of the corrected road node on the bridge
  pt_junction <- apply(tb_node_bridge[,c("X","Y")], 2, mean) |>
    sf::st_point() |>
    sf::st_sfc(crs = sf::st_crs(roads))

  # Construct bridge with coordinates in the right order
  bridge <- rbind(coord_1[-nrow(coord_1),], sf::st_coordinates(pt_junction), coord_2[-1,]) |>
    sf::st_linestring() |>
    sf::st_sfc(crs = sf::st_crs(roads))

  return(list(bridge = bridge,
              junction = pt_junction,
              reversed = c(rev_1, rev_2)))
}