Refactor rnet_breakup_vertices (#451)

* First ideas for #416

* fix problem in the interaction between route() and data.table created after setting  .datatable.aware = TRUE

* Second round

* Improve rnet_breakup

* fix problem with CRS in rnet_breakup

DESCRIPTION

@@ -51,7 +51,9 @@ Imports:
     rlang (>= 0.2.2),
     lwgeom (>= 0.1.4),
     sf (>= 0.6.3),
-    magrittr
+    magrittr,
+    sfheaders,
+    data.table
 LinkingTo: 
     RcppArmadillo (>= 0.9.100.5.0),
     Rcpp (>= 0.12.18)
@@ -70,11 +72,9 @@ Suggests:
     bench,
     openxlsx (>= 4.1.0),
     osrm,
-    data.table,
     geodist,
     mapsapi,
-    s2,
-    sfheaders
+    s2
 VignetteBuilder: knitr
 URL: https://github.com/ropensci/stplanr, https://docs.ropensci.org/stplanr/
 SystemRequirements: GNU make

R/rnet-clean.R

@@ -1,36 +1,39 @@
-#' Break up an `sf` object with LINESTRING geometry.
+# See https://rdatatable.gitlab.io/data.table/articles/datatable-importing.html#data-table-in-imports-but-nothing-imported-1
+.datatable.aware = TRUE
+utils::globalVariables(c("linestring_id", "new_linestring_id"))
+
+#' Break up an sf object with LINESTRING geometry.
 #'
-#' This function breaks up a single LINESTRING geometry into multiple
-#' LINESTRING(s) for preserving routability of an `sfNetwork` object created
-#' by [SpatialLinesNetwork()] function with Open Street Map data. See details
-#' and [stplanr/issues/282](https://github.com/ropensci/stplanr/issues/282).
+#' This function breaks up a LINESTRING geometry into multiple LINESTRING(s). It
+#' is used mainly for preserving routability of an `sfNetwork` object that is
+#' created using Open Street Map data. See details,
+#' [stplanr/issues/282](https://github.com/ropensci/stplanr/issues/282), and
+#' [stplanr/issues/416](https://github.com/ropensci/stplanr/issues/416).
 #'
-#' A LINESTRING geometry is broken-up when one of the following conditions is
-#' met:
-#' 1. two or more LINESTRINGS share a POINT that lies in the union of their
-#' boundaries (see the rnet_roundabout example);
+#' A LINESTRING geometry is broken-up when one of the two following conditions
+#' are met:
+#' 1. two or more LINESTRINGS share a POINT which is a boundary point for some
+#' LINESTRING(s), but not all of them (see the rnet_roundabout example);
 #' 2. two or more LINESTRINGS share a POINT which is not in the boundary of any
 #' LINESTRING (see the rnet_cycleway_intersection example).
 #'
 #' The problem with the first example is that, according to algorithm behind
-#' [SpatialLinesNetwork()], two LINESTRINGS are connected if and only
-#' if they share at least one point in their boundaries. The roads and the
-#' roundabout are clearly connected in the "real" world but the corresponding
-#' LINESTRING objects do not share any boundary point. In fact, by Open Street
+#' [SpatialLinesNetwork()], two LINESTRINGS are connected if and only if they
+#' share at least one point in their boundaries. The roads and the roundabout
+#' are clearly connected in the "real" world but the corresponding LINESTRING
+#' objects do not share two distinct boundary points. In fact, by Open Street
 #' Map standards, a roundabout is represented as a closed and circular
-#' LINESTRING and this implies that the roundabout is not connected to the other
-#' roads according to [SpatialLinesNetwork()] definition. By the same reasoning,
-#' the roads in the second example are clearly connected in the "real" world,
-#' but they do not share any point in their boundaries. This function is used to
-#' solve this type of problem.
+#' LINESTRING, and this implies that the roundabout is not connected to the
+#' other roads according to [SpatialLinesNetwork()] definition. By the same
+#' reasoning, the roads in the second example are clearly connected in the
+#' "real" world, but they do not share any point in their boundaries. This
+#' function is used to solve this type of problem.
 #'
-#' @param rnet An sf object with LINESTRING geometry representing a route
+#' @param rnet An sf or sfc object with LINESTRING geometry representing a route
 #'   network.
-#' @param breakup_internal_vertex_matches Boolean. Should breaks be made at
-#'   shared internal points? `TRUE` by default. Internal points are points that
-#'   do not lie in the boundary of the LINESTRING.
-#' @return An sf object with LINESTRING geometry created after breaking up the
-#'   input object.
+#' @param verbose Boolean. If TRUE, the function prints additional messages.
+#' @return An sf or sfc object with LINESTRING geometry created after breaking
+#'   up the input object.
 #' @family rnet
 #' @export
 #'
@@ -40,13 +43,13 @@
 #' par(mar = rep(0, 4))
 #'
 #' # Check the geometry of the roundabout example. The dots represent the
-#' # boundary points of the LINESTRINGS. The "isolated" red point in the top-left
-#' # is the boundary point of the roundabout, and it is not shared with any
-#' # other street.
+#' # boundary points of the LINESTRINGS. The "isolated" red point in the
+#' # top-left is the boundary point of the roundabout, and it is not shared
+#' # with any other street.
 #' plot(st_geometry(rnet_roundabout), lwd = 2, col = rainbow(nrow(rnet_roundabout)))
 #' boundary_points <- st_geometry(line2points(rnet_roundabout))
 #' points_cols <- rep(rainbow(nrow(rnet_roundabout)), each = 2)
-#' plot(boundary_points, pch = 16, add = TRUE, col = points_cols)
+#' plot(boundary_points, pch = 16, add = TRUE, col = points_cols, cex = 2)
 #'
 #' # Clean the roundabout example.
 #' rnet_roundabout_clean <- rnet_breakup_vertices(rnet_roundabout)
@@ -56,14 +59,14 @@
 #' plot(boundary_points, pch = 16, add = TRUE, col = points_cols)
 #' # The roundabout is now routable since it was divided into multiple pieces
 #' # (one for each colour), which, according to SpatialLinesNetwork() function,
-#' # are connected to the other streets.
+#' # are connected.
 #'
 #' # Check the geometry of the overpasses example. This example is used to test
 #' # that this function does not create any spurious intersection.
 #' plot(st_geometry(rnet_overpass), lwd = 2, col = rainbow(nrow(rnet_overpass)))
 #' boundary_points <- st_geometry(line2points(rnet_overpass))
 #' points_cols <- rep(rainbow(nrow(rnet_overpass)), each = 2)
-#' plot(boundary_points, pch = 16, add = TRUE, col = points_cols)
+#' plot(boundary_points, pch = 16, add = TRUE, col = points_cols, cex = 2)
 #' # At the moment the network is not routable since one of the underpasses is
 #' # not connected to the other streets.
 #'
@@ -80,61 +83,145 @@
 #' # Check the geometry of the cycleway_intersection example. The black dots
 #' # represent the boundary points and we can see that the two roads are not
 #' # connected according to SpatialLinesNetwork() function.
-#' plot(rnet_cycleway_intersection$geometry,
+#' plot(
+#'   rnet_cycleway_intersection$geometry,
 #'   lwd = 2,
-#'   col = rainbow(nrow(rnet_cycleway_intersection))
+#'   col = rainbow(nrow(rnet_cycleway_intersection)),
+#'   cex = 2
 #' )
 #' plot(st_geometry(line2points(rnet_cycleway_intersection)), pch = 16, add = TRUE)
 #' # Check interactively
 #' # mapview::mapview(rnet_overpass)
 #'
 #' # Clean the rnet object and plot the result.
 #' rnet_cycleway_intersection_clean <- rnet_breakup_vertices(rnet_cycleway_intersection)
-#' plot(rnet_cycleway_intersection_clean$geometry,
-#'   lwd = 2, col = rainbow(nrow(rnet_cycleway_intersection_clean))
+#' plot(
+#'   rnet_cycleway_intersection_clean$geometry,
+#'   lwd = 2,
+#'   col = rainbow(nrow(rnet_cycleway_intersection_clean)),
+#'   cex = 2
 #' )
 #' plot(st_geometry(line2points(rnet_cycleway_intersection_clean)), pch = 16, add = TRUE)
 #'
 #' par(def_par)
-rnet_breakup_vertices <- function(rnet, breakup_internal_vertex_matches = TRUE) {
-  rnet_nodes <- sf::st_geometry(line2points(rnet))
-  rnet_internal_vertexes <- sf::st_geometry(line2vertices(rnet))
-
-  # For the first part of the procedure I don't need duplicated nodes or
-  # duplicated vertexes so I can extract their unique values
-  unique_rnet_nodes <- do.call("c", unique(rnet_nodes))
-  unique_rnet_internal_vertexes <- do.call("c", unique(rnet_internal_vertexes))
-
-  # Intersection between nodes and internal vertexes
-  # The following code is the same as
-  # intersection_point <- sf::st_intersection(unique_rnet_nodes, unique_rnet_internal_vertexes)
-  # but faster since we are dealing only with points
-
-  rbind_nodes_internal_vertexes <- rbind(unique_rnet_nodes, unique_rnet_internal_vertexes)
-  index_intersection_points <- duplicated(rbind_nodes_internal_vertexes)
-
-  if (any(index_intersection_points)) {
-    intersection_points <- sf::st_as_sf(
-      data.frame(rbind_nodes_internal_vertexes[index_intersection_points, , drop = FALSE]),
-      coords = c("x_coords", "y_coords"),
-      crs = sf::st_crs(rnet)
+
+rnet_breakup_vertices <- function(rnet, verbose = FALSE) {
+  # A few safety checks
+  if (!inherits(rnet, c("sf", "sfc"))) {
+    stop(
+      "Sorry, at the moment this function works only with sf and sfc objects",
+      call. = FALSE
+    )
+  }
+  if (!all(sf::st_is(rnet, "LINESTRING"))) {
+    stop(
+      "The input object must have LINESTRING geometry. See also ?sf::st_cast.",
+      call. = FALSE
     )
+  }
+
+  # Check if we need to rebuild the tbl_df structure at the end
+  rebuild_tbl <- FALSE
+  if (inherits(rnet, "sf") && inherits(rnet, "tbl_df")) {
+    rebuild_tbl <- TRUE
+  }
 
-    message("Splitting rnet object at the shared boundary points.")
-    rnet_breakup_collection <- lwgeom::st_split(rnet, intersection_points$geometry)
-    rnet_clean <- sf::st_collection_extract(rnet_breakup_collection, "LINESTRING")
-  } else {
-    rnet_clean <- rnet
+  # Step 1 - We need to split all LINESTRING(s) using the points that are both
+  # duplicated points and not in the boundary
+
+  # 1a - Extract all points and convert to data.table format
+  rnet_points <- data.table::data.table(
+    sfheaders::sfc_to_df(sf::st_geometry(rnet))
+  )
+  # The existing_dimensions vector is used to perform the following operations using x,y; x,y,z and so on.
+  existing_dimensions <- intersect(c("x", "y", "z", "m"), colnames(rnet_points))
+
+  # 1b - Extract id of boundary points
+  id_boundary_points <- rnet_points[
+    !duplicated(linestring_id) | !duplicated(linestring_id, fromLast = TRUE),
+    which = TRUE
+  ]
+  if (verbose) {
+    message("Extracted the ID(s) of the boundary points")
   }
 
-  # Split again at the duplicated internal vertexes
-  rnet_internal_vertexes_duplicated <- rnet_internal_vertexes[duplicated(rnet_internal_vertexes)]
+  # 1c - Subset all duplicated points
+  id_duplicated_points <- rnet_points[
+    duplicated(rnet_points, by = existing_dimensions) |
+    duplicated(rnet_points, by = existing_dimensions, fromLast = TRUE),
+    which = TRUE
+  ]
+  if (verbose) {
+    message("Extracted the ID(s) of the duplicated points")
+  }
+
+  # 1d - Anti-join to find the set difference between duplicated_points and
+  # boundary_points
+  id_split_points <- setdiff(id_duplicated_points, id_boundary_points)
+
+  if (length(id_split_points) == 0L) {
+    message(
+      "The input data doesn't have any duplicated point which is also a boundary point.",
+      "Return the input data."
+    )
+    return(rnet)
+  }
+
+  # Now I need to duplicate the coordinates of the internal points that are also
+  # boundary points (since they will become boundary points for the new
+  # LINESTRING(s)). I use sort since I want to preserve the order within the
+  # points.
+  rnet_points <- rnet_points[
+    sort(c(seq_len(nrow(rnet_points)), id_split_points))
+  ]
+
+  # I will build the new LINESTRING using sfheaders::sfc_linestring so I need to
+  # build an ID for each new LINESTRING. The ID for the new LINESTRING(s) is
+  # created starting from the linestring_id column in rnet_points and
+  # incrementing the id by 1 at each break point.
+  break_id <- double(nrow(rnet_points))
+  # The ID is shifted by 1 at each break point
+  break_id[id_split_points + 1:length(id_split_points)] <- 1L
+  break_id <- cumsum(break_id)
+  rnet_points[["new_linestring_id"]] <- rnet_points[["linestring_id"]] + break_id
+
+  # Now I can create the new LINESTRING sfc
+  keep_columns <- c(existing_dimensions, "new_linestring_id")
+  new_linestring <- sfheaders::sfc_linestring(
+    rnet_points[, keep_columns, with = FALSE],
+    linestring_id = "new_linestring_id"
+  )
+
+  # Add CRS and precision
+  new_linestring_sfc <- sf::st_sfc(
+    new_linestring,
+    crs = sf::st_crs(rnet),
+    precision = sf::st_precision(rnet)
+  )
+
+  # If the input is sfc then I just need to return new_linestring_sfc, otherwise
+  # I have to rebuild the sf structure.
+  if (!inherits(rnet, "sf")) {
+    return(new_linestring_sfc)
+  }
+
+  # Determine the ID(s) of the old LINESTRING(s) (i.e. I need to create an sf
+  # object with the new geometry column and the old fields)
+  old_rnet_id <- rnet_points[
+    j = list(linestring_id = unique(linestring_id)),
+    by = new_linestring_id
+  ][["linestring_id"]]
+
+  rnet <- sf::st_sf(
+    sf::st_drop_geometry(rnet)[old_rnet_id, , drop = FALSE],
+    geometry = new_linestring_sfc,
+    agr = sf::st_agr(rnet)
+  )
 
-  if (length(rnet_internal_vertexes_duplicated) > 0 & breakup_internal_vertex_matches) {
-    message("Splitting rnet object at the shared internal points.")
-    rnet_breakup_collection <- lwgeom::st_split(rnet_clean, rnet_internal_vertexes_duplicated)
-    rnet_clean <- sf::st_collection_extract(rnet_breakup_collection, "LINESTRING")
+  # 5 - Maybe we need to rebuild the tbl_df structure
+  if (rebuild_tbl) {
+    rnet <- sf::st_as_sf(dplyr::as_tibble(rnet))
   }
 
-  rnet_clean
+  rnet
 }

R/route.R

@@ -110,7 +110,7 @@ route.sf <- function(from = NULL, to = NULL, l = NULL,
     # browser()
     # warning("data.table used to create the sf object, bounding box may be incorrect.")
     out_dt <- data.table::rbindlist(list_out[list_elements_sf])
-    out_dtsf <- sf::st_sf(out_dt[, !names(out_dt) %in% "geometry"], geometry = out_dt$geometry)
+    out_dtsf <- sf::st_sf(out_dt[, !"geometry"], geometry = out_dt$geometry)
     # attributes(out_dtsf$geometry)
     # identical(sf::st_bbox(out_dtsf), sf::st_bbox(out_sf)) # FALSE
     attr(out_dtsf$geometry, "bbox") = sfheaders::sf_bbox(out_dtsf)