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:/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:/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:/ropensci/stplanr/issues/282), and
+#' [stplanr/issues/416](https:/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)