Merge branch 'main' into 11-tests-cf

DESCRIPTION

@@ -13,25 +13,27 @@ Description: CRiSp (City River Spaces) provides tools to automate the
 License: Apache License (>= 2)
 URL: https://cityriverspaces.github.io/CRiSp/
 BugReports: https://github.com/CityRiverSpaces/crisp/issues
+Depends: 
+    R (>= 2.10)
 Imports: 
     dplyr,
     lwgeom,
     osmdata,
+    rlang,
     sf,
     sfnetworks,
     stringr,
-    tidygraph,
-    rlang
+    tidygraph
 Suggests: 
+    ggplot2,
+    gridExtra,
     knitr,
     rmarkdown,
     testthat (>= 3.0.0)
 VignetteBuilder: 
     knitr
 Config/testthat/edition: 3
 Encoding: UTF-8
+LazyData: true
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.3.2
-Depends: 
-    R (>= 2.10)
-LazyData: true

NAMESPACE

@@ -6,12 +6,14 @@ export(clean_network)
 export(create_network)
 export(define_aoi)
 export(delineate_corridor)
+export(flatten_network)
 export(get_corridor_edge)
-export(get_latlon)
+export(get_geom_latlon)
 export(get_osm_city_boundary)
 export(get_osmdata)
 export(get_osmdata_river_corridor)
 export(get_target_points)
+export(get_utm_zone_epsg)
 export(get_vertices)
 export(merge_streets)
 export(not_intersects)

R/delineate-corridor.R

@@ -55,6 +55,9 @@ trim_network <- function(net, area, river_corridor) {
 
 #' Simplify a street network by removing multiple edges and loops.
 #'
+#' Simplify the graph, removing loops and double-edge connections following
+#' [this approach](https://luukvdmeer.github.io/sfnetworks/articles/sfn02_preprocess_clean.html#simplify-network).
+#'
 #' @param net A network object
 #'
 #' @return A simplifed network object
@@ -71,15 +74,21 @@ simplify_network <- function(net) {
 
 #' Clean a street network by subdividing edges and removing pseudo-nodes.
 #'
+#' Subdivide edges by [adding missing nodes](https://luukvdmeer.github.io/sfnetworks/articles/sfn02_preprocess_clean.html#subdivide-edges),
+#' remove [pseudo-nodes](https://luukvdmeer.github.io/sfnetworks/articles/sfn02_preprocess_clean.html#smooth-pseudo-nodes),
+#' and keep only the main connected component of the network.
+#'
 #' @param net A network object
 #'
 #' @return A cleaned network object
 #' @export
 clean_network <- function(net) {
   net |>
+    tidygraph::convert(sfnetworks::to_spatial_subdivision, .clean = TRUE) |>
+    tidygraph::convert(sfnetworks::to_spatial_smooth, .clean = TRUE) |>
     CRiSp::simplify_network() |>
-    tidygraph::convert(sfnetworks::to_spatial_subdivision) |>
-    tidygraph::convert(sfnetworks::to_spatial_smooth)
+    tidygraph::activate("nodes") |>
+    dplyr::filter(tidygraph::group_components() == 1)
 }
 
 #' Determine the end vertices of the initial river corridor.

R/preprocess.R

@@ -31,3 +31,127 @@ create_network <- function(data, crs = NULL) {
 
   net
 }
+
+#' Flatten a network by adding points at apparent intersections
+#'
+#' This function adds nodes at intersections between network edges that
+#' cross each other but do not have vertices at the intersection.
+#'
+#' @param network A network object
+#'
+#' @return A network object with additional points at intersections
+#' @export
+flatten_network <- function(network) {
+  # Determine intersection points between crossing edges
+  edges_cross <- get_crossing_edges(network)
+  pts_intersect <- get_intersection_points(edges_cross)
+
+  # Convert edge table to data.frame and add info on boundary points
+  edge_pts <- sfheaders::sf_to_df(edges_cross)
+  edge_idxs <- edge_pts$linestring_id
+  edge_pts$is_startpoint <- !duplicated(edge_idxs)
+  edge_pts$is_endpoint <- !duplicated(edge_idxs, fromLast = TRUE)
+
+  # Loop over all points, add them to the edge table
+  for (i in seq_len(nrow(pts_intersect))) {
+    point <- pts_intersect$geometry[[i]]
+    intersecting_edges <- unique(unlist(pts_intersect$origins[i]))
+    for (edge_id in intersecting_edges){
+      edge_pts <- insert_intersection(edge_pts, point, edge_id)
+    }
+  }
+
+  # Convert back edge table to sfc object
+  edges_cross_new <- sfheaders::sfc_linestring(edge_pts, linestring_id = "id",
+                                               x = "x", y = "y")
+  sf::st_crs(edges_cross_new) <- sf::st_crs(edges_cross)
+
+  # Update the network with the new edge geometries
+  nodes <- network |> sf::st_as_sf("nodes")
+  edges <- network |> sf::st_as_sf("edges")
+  edges[edges_cross$id, ] <- edges[edges_cross$id, ] |>
+    sf::st_set_geometry(edges_cross_new)
+  network_new <- sfnetworks::sfnetwork(
+    nodes = nodes,
+    edges = edges,
+    directed = FALSE,
+    force = TRUE,  # skip checks
+  )
+  network_new
+}
+
+get_crossing_edges <- function(network) {
+  network |>
+    tidygraph::activate("edges") |>
+    # Add ID to ease replacement later on
+    dplyr::mutate(id = seq_len(dplyr::n())) |>
+    dplyr::filter(sfnetworks::edge_crosses(tidygraph::.E())) |>
+    sf::st_as_sf("edges")
+}
+
+get_intersection_points <- function(edges) {
+  pts_intersect <- edges |>
+    sf::st_intersection() |>
+    # Cast multipoint intersections into points
+    sf::st_collection_extract("POINT") |>
+    sfheaders::sf_cast(to = "POINT")
+
+  pts_intersect_agg <- aggregate(
+    pts_intersect,
+    by = sf::st_geometry(pts_intersect),
+    FUN = unique,
+    drop = TRUE
+  )
+
+  pts_intersect_unique <- pts_intersect_agg |> dplyr::distinct()
+  pts_intersect_unique
+}
+
+distance <- function(x1, y1, x2, y2) {
+  sqrt((x2 - x1) ^ 2 + (y2 - y1) ^ 2)
+}
+
+insert_intersection <- function(edge_pts, point, line_id) {
+  line_pts <- subset(edge_pts, linestring_id == line_id)
+  pt_x <- point[[1]]
+  pt_y <- point[[2]]
+  is_point_in_line <- nrow(subset(line_pts, x == pt_x &
+                                    y == pt_y)) >= 1
+  if (!is_point_in_line) {
+    startpoint <- subset(line_pts, is_startpoint == TRUE)
+    kk <- as.numeric(rownames(startpoint))
+    w_break <- FALSE
+    while (!w_break) {
+      # Consider the line segments a - b.
+      # x is a valid intersection if the following condition is true:
+      # distance(a, b) == distance(a, x) + distance(x, b)
+      pt_a_x <- edge_pts[kk, ]$x
+      pt_a_y <- edge_pts[kk, ]$y
+      pt_b_x <- edge_pts[kk + 1, ]$x
+      pt_b_y <- edge_pts[kk + 1, ]$y
+      d_ab <- distance(pt_a_x, pt_a_y, pt_b_x, pt_b_y)
+      d_ax <- distance(pt_a_x, pt_a_y, pt_x, pt_y)
+      d_bx <- distance(pt_b_x, pt_b_y, pt_x, pt_y)
+      is_intersection <- dplyr::near(d_ab, d_ax + d_bx, tol = 1.e-3)
+      if (is_intersection) {
+        insertion <- tibble::tibble_row(
+          sfg_id = line_id,
+          linestring_id = line_id,
+          x = pt_x,
+          y = pt_y,
+          is_startpoint = FALSE,
+          is_endpoint = FALSE
+        )
+        edge_pts <- tibble::add_row(edge_pts, insertion, .after = kk)
+        w_break <- TRUE
+      } else {
+        if (edge_pts[kk + 1, ]$is_endpoint) {
+          warning("point is not added to the edge df.")
+          w_break <- TRUE
+        }
+      }
+      kk <- kk + 1
+    }
+  }
+  edge_pts
+}

R/utils.R

@@ -6,7 +6,7 @@
 #'
 #' @return An sf object with latitude and longitude
 #' @export
-get_latlon <- function(x) {
+get_geom_latlon <- function(x) {
   sf::st_transform(x, 4326) |>
     sf::st_geometry()
 }
@@ -34,3 +34,29 @@ calc_weights <- function(net) {
     sfnetworks::activate("edges") |>
     dplyr::mutate(weight = sfnetworks::edge_length())
 }
+
+#' Get UTM zone from longitude
+#'
+#' @param x An sf object
+#'
+#' @return The UTM zone
+#' @export
+get_utm_zone_epsg <- function(x) {
+  if (!"sf" %in% class(x)) {
+    stop("x must be an sf object")
+  }
+
+  coords <- x |>
+    sf::st_bbox() |>
+    sf::st_as_sfc() |>
+    sf::st_centroid() |>
+    sf::st_transform(4326) |>
+    sf::st_coordinates()
+
+  if (coords[2] >= 0L) {
+    base <- 32600
+  } else {
+    base <- 32700
+  }
+  base + floor((coords[1] + 180) / 6) + 1
+}

_pkgdown.yml

@@ -9,6 +9,8 @@ articles:
   navbar: ~
   contents:
   - method
+  - getting-osm-data
+  - network-preparation
   - corridor-delineation
   - corridor-segmentation
   - riverspace-delineation
@@ -21,15 +23,23 @@ navbar:
       menu:
         - text: 1. The method
           href: articles/method.html
-        - text: 2. Corridor delineation
+        - text: -------
+        - text: Preparation
+        - text: 2. Getting OSM data for delineation
+          href: articles/getting-osm-data.html
+        - text: 3. Preparing the network for delineation
+          href: articles/network-preparation.html
+        - text: -------
+        - text: Delineation
+        - text: 4. Corridor delineation
           href: articles/corridor-delineation.html
-        - text: 3. Corridor segmentation
+        - text: 5. Corridor segmentation
           href: articles/corridor-segmentation.html
-        - text: 4. Riverspace delineation
+        - text: 6. Riverspace delineation
           href: articles/riverspace-delineation.html
         - text: -------
         - text: Use cases
-        - text: 5. Multiple delineations
+        - text: 7. Multiple delineations
           href: articles/multiple-cities.html
-        - text: 6. Study area around a POI
+        - text: 8. Study area around a POI
           href: articles/poi-study-area.html

data-raw/bucharest.R

@@ -10,15 +10,20 @@ bbox_buffer <- 2000
 
 city_boundary <- get_osm_city_boundary(city_name) |>
   st_transform(epsg_code)
+st_crs(city_boundary)$wkt <- gsub("°|º", "\\\u00b0",
+                                  st_crs(city_boundary)$wkt)
 
 bb <- getbb(city_name)
 aoi <- define_aoi(bb, epsg_code, bbox_buffer)
+bbox_expanded <- aoi |> st_transform(4326) |> st_bbox()
 
 river_centerline <- osmdata_as_sf("waterway", "river", bb)$osm_multilines |>
   filter(name == river_name) |>
   st_transform(epsg_code) |>
   st_geometry() |>
   st_intersection(st_buffer(aoi, bbox_buffer))
+st_crs(river_centerline)$wkt <- gsub("°|º", "\\\u00b0",
+                                     st_crs(river_centerline)$wkt)
 
 river_surface <- osmdata_as_sf("natural", "water", bb)
 river_surface <- river_surface$osm_multipolygons |>
@@ -27,18 +32,31 @@ river_surface <- river_surface$osm_multipolygons |>
   st_filter(river_centerline, .predicate = st_intersects) |>
   st_geometry() |>
   st_union()
+st_crs(river_surface)$wkt <- gsub("°|º", "\\\u00b0",
+                                  st_crs(river_surface)$wkt)
 
 highway_values <- c("motorway", "primary", "secondary", "tertiary")
 streets <- osmdata_as_sf("highway", highway_values, bb)
 streets <- merge_streets(streets) |>
-  select("highway")
+  select("highway") |>
+  st_transform(epsg_code)
+st_crs(streets)$wkt <- gsub("°|º", "\\\u00b0", st_crs(streets)$wkt)
+
+railways <- osmdata_as_sf("railway", "rail", bbox_expanded)
+railways_lines <- railways$osm_lines |>
+  select("railway") |>  # only keep "railway" column
+  rename(type = `railway`) |>  # rename it to "type"
+  st_transform(epsg_code)
+st_crs(railways_lines)$wkt <- gsub("°|º", "\\\u00b0",
+                                   st_crs(railways_lines)$wkt)
 
 bucharest <- list(
   bb = bb,
   boundary = city_boundary,
   river_centerline = river_centerline,
   river_surface = river_surface,
-  streets = streets
+  streets = streets,
+  railways_lines = railways_lines
 )
 
 usethis::use_data(bucharest, overwrite = TRUE)