From e9f54d607b0b4bf859431e13822d1b9affa01906 Mon Sep 17 00:00:00 2001 From: Eduardo Leoni Date: Fri, 4 Apr 2025 10:13:51 -0400 Subject: [PATCH 1/4] add_nodes_to_graph_by_edge proposal --- R/match-points-by-edge.R | 325 ++++++++++++++++++++ tests/testthat/test-add-nodes-by-edge.R | 378 ++++++++++++++++++++++++ 2 files changed, 703 insertions(+) create mode 100644 R/match-points-by-edge.R create mode 100644 tests/testthat/test-add-nodes-by-edge.R diff --git a/R/match-points-by-edge.R b/R/match-points-by-edge.R new file mode 100644 index 00000000..86b5a02b --- /dev/null +++ b/R/match-points-by-edge.R @@ -0,0 +1,325 @@ +#' Add nodes to a graph by matching them to the nearest edges, processing +#' multiple points on the same edge efficiently. +#' +#' This function extends \link{add_nodes_to_graph} by handling cases where multiple +#' points match to the same edge more efficiently. It splits edges at points of +#' perpendicular intersection and optionally creates edges connecting to the original +#' points. +#' +#' @inheritParams add_nodes_to_graph +#' @param xy_id Optional vector of IDs to assign to the new nodes. If NULL, IDs +#' will be generated automatically. +#' @param wt_profile Character string specifying the weight profile to use for +#' edges connecting to original points (e.g., "foot", "bike", "car"). +#' @param wt_profile_file Character string specifying the file path to a custom +#' weight profile CSV file. Used only if `wt_profile` is NULL. +#' @param highway Character string specifying the highway type for edges connecting +#' to original points. Must be a type defined in the weight profile. +#' @param max_distance Numeric value specifying the maximum allowed distance for edges +#' connecting to original points. Points beyond this distance won't be connected. +#' @param replace_component Logical value specifying whether to recalculate graph +#' component IDs after adding the new nodes and edges. +#' @return Modified version of graph with nodes added at specified locations, +#' original edges split at intersection points, and (unless `intersections_only=TRUE`) +#' additional edges connecting to the original points. +#' @export +#' +add_nodes_to_graph_by_edge <- function (graph, + xy, + dist_tol = 1e-6, + intersections_only = FALSE, + xy_id = NULL, + wt_profile = NULL, + wt_profile_file = NULL, + highway = NULL, + max_distance = Inf, + replace_component = TRUE) { + + genhash <- function (len = 10) { + paste0 (sample (c (0:9, letters, LETTERS), size = len), collapse = "") + } + genhashv <- Vectorize (genhash) + + xy <- pre_process_xy (xy) + + gr_cols <- dodgr_graph_cols (graph) + gr_cols <- unlist (gr_cols [which (!is.na (gr_cols))]) + graph_std <- graph [, gr_cols] # standardise column names + names (graph_std) <- names (gr_cols) + + # Match points to the standardized graph + pts <- match_pts_to_graph (graph_std, xy, distances = TRUE) + projids <- unique (pts [, c ("x", "y")]) + projids$proj_id <- paste0 ("proj_", genhash (), "_", seq_len (nrow (projids))) + pts <- merge (pts, projids) + + if (is.null (xy_id)) { + pts$xy_id <- genhashv (rep (10L, nrow (pts))) + } else { + pts$xy_id <- xy_id + } + + # Return original graph if no points match + if (nrow (pts) == 0) { + return (graph) + } + + graph$tmp_graph_index <- 1:nrow (graph) + + # Add original coordinates to pts + pts$x0 <- xy [, 1] + pts$y0 <- xy [, 2] + pts$pts_index <- seq_len (nrow (pts)) + pts$from <- graph_std$from [pts$index] + pts$to <- graph_std$to [pts$index] + + + # Create a lookup table from the graph for efficient matching + graph_lookup <- data.frame ( + from = graph_std$from, + to = graph_std$to, + index = seq_len (nrow (graph_std)) + ) + + # Create bidirectional version with swapped from/to + pts_bi <- pts [, setdiff (names (pts), "index")] # Remove index column + + # Find matching indices using merge with swapped column order + pts_bi <- merge ( + graph_lookup, + pts_bi, + by.y = c ("to", "from"), # Swap the order in the merge operation + by.x = c ("from", "to") + ) + + # Combine with original points + pts <- unique (rbind (pts, pts_bi [names (pts)])) + + # Extract edges that need to be split + edges_to_split <- graph_std [pts$index, ] + + # Add index for tracking + edges_to_split$n <- pts$pts_index + + graph_to_add <- graph [pts$index, ] + + # Remove edges to be split from the graph + graph_std <- graph_std [-pts$index, ] + graph <- graph [-pts$index, ] + + # Group edges by edge_id + unique_edges <- unique (edges_to_split$edge_id) + all_edges_split <- list () + + # Process each unique edge + for (edge_id in unique_edges) { + # Get all instances of this edge + current_edge <- edges_to_split [edges_to_split$edge_id == edge_id, ] + current_edge_1 <- current_edge [1, ] + + # Get all points that match this edge instance + edge_pts <- pts [pts$pts_index %in% current_edge$n, ] + edge_pts <- edge_pts [!duplicated (edge_pts$proj_id), ] + + # stop if no points + if (nrow (edge_pts) == 0) { + stop ("No points found.") + } + + # Store the original ratios - these will be used directly + orig_d_weighted <- current_edge_1$d_weighted + orig_d <- current_edge_1$d + orig_ratio <- orig_d_weighted / orig_d + orig_time_ratio <- current_edge_1$time_weighted / current_edge_1$time + + # Sort points along the edge + # Calculate distance from start of edge to each projection point + start_point <- c (current_edge$xfr [1], current_edge$yfr [1]) + end_point <- c (current_edge$xto [1], current_edge$yto [1]) + + # Calculate vector from start to end + edge_vector <- end_point - start_point + edge_length <- sqrt (sum (edge_vector^2)) + + # Calculate projection of each point onto the edge + proj_distances <- numeric (nrow (edge_pts)) + for (p in seq_len (nrow (edge_pts))) { + point_vector <- c (edge_pts$x [p], edge_pts$y [p]) - start_point + # Projection distance along the edge + proj_distances [p] <- sum (point_vector * edge_vector) / edge_length + } + + # Sort points by projection distance + sorted_indices <- order (proj_distances) + edge_pts <- edge_pts [sorted_indices, ] + proj_distances <- proj_distances [sorted_indices] + proj_ok <- which (proj_distances > dist_tol) + if (length (proj_ok) == 0) { + proj_ok <- 1 + } + edge_pts <- edge_pts [proj_ok, ] + + # Process all points normally without special handling for close points + # Create segments for each point + all_segments <- NULL + + # If not intersections_only, add connections to original points + if (!intersections_only) { + for (p in seq_len (nrow (edge_pts))) { + # Create edges to original point (bidirectional) + new_edges <- rbind (current_edge_1, current_edge_1) + + # Set up connection to original point + new_edges$from [1] <- new_edges$to [2] <- edge_pts$xy_id [p] + new_edges$to [1] <- new_edges$from [2] <- edge_pts$proj_id [p] + new_edges$xfr [1] <- new_edges$xto [2] <- edge_pts$x0 [p] + new_edges$yfr [1] <- new_edges$yto [2] <- edge_pts$y0 [p] + new_edges$xto [1] <- new_edges$xfr [2] <- edge_pts$x [p] + new_edges$yto [1] <- new_edges$yfr [2] <- edge_pts$y [p] + + # Calculate distance using geodesic distance + d_i <- geodist::geodist ( + data.frame ( + x = c (new_edges$xfr [1], new_edges$xto [1]), + y = c (new_edges$yfr [1], new_edges$yto [1]) + ), + measure = "geodesic" + ) [1, 2] + + # Skip if distance smaller than dist_tol + if (d_i < dist_tol) { + next + } + + # Skip if distance exceeds maximum allowed + if (d_i > max_distance) { + next + } + + # Apply custom weight profile if provided + if (!is.null (wt_profile) || !is.null (wt_profile_file)) { + # Get weight profile + wp <- dodgr:::get_profile (wt_profile = wt_profile, file = wt_profile_file) + way_wt <- wp$value [wp$way == highway] + + if (length (way_wt) == 0) { + stop ("Highway type '", highway, "' not found in weight profile") + } + + # Calculate weights using the profile + new_edges$d <- d_i + new_edges$d_weighted <- d_i / way_wt + new_edges$highway <- highway + # Apply additional weighting functions + new_edges <- dodgr:::set_maxspeed (new_edges, wt_profile, wt_profile_file) |> + dodgr:::weight_by_num_lanes (wt_profile) |> + dodgr:::calc_edge_time (wt_profile) + + } else { + # Use original edge's weight ratios + new_edges$d <- d_i + new_edges$d_weighted <- d_i * orig_ratio + new_edges$time <- d_i * (current_edge_1$time / current_edge_1$d) + new_edges$time_weighted <- d_i * (current_edge_1$time / current_edge_1$d) * orig_time_ratio + new_edges$highway <- unique (graph_to_add [graph_to_add$tmp_graph_index == edge_pts$index [p], "highway"]) + } + + # Generate unique edge IDs + new_edges$edge_id <- c ( + paste0 ("to_orig_", p, "_", genhash ()), + paste0 ("from_orig_", p, "_", genhash ()) + ) + + # Add to all segments + all_segments <- rbind (all_segments, new_edges) + } + } + + # Now split the edge into segments + n_segments <- nrow (edge_pts) + 1 + + # Create a list to store all segments + segments <- vector ("list", n_segments) + + # Create first segment + segments [[1]] <- current_edge_1 + segments [[1]]$to <- edge_pts$proj_id [1] + segments [[1]]$xto <- edge_pts$x [1] + segments [[1]]$yto <- edge_pts$y [1] + + # Create middle segments + if (n_segments > 2) { + for (s in seq_len (n_segments - 2)) { + segments [[s + 1]] <- current_edge_1 + segments [[s + 1]]$from <- segments [[s]]$to + segments [[s + 1]]$to <- edge_pts$proj_id [s] + segments [[s + 1]]$xfr <- edge_pts$x [s] + segments [[s + 1]]$yfr <- edge_pts$y [s] + segments [[s + 1]]$xto <- edge_pts$x [s + 1] + segments [[s + 1]]$yto <- edge_pts$y [s + 1] + } + } + + # Create last segment + segments [[n_segments]] <- current_edge_1 + segments [[n_segments]]$from <- if (n_segments > 1) { + segments [[n_segments - 1]]$to + } else { + segments [[1]]$to + } + segments [[n_segments]]$xfr <- edge_pts$x [n_segments - 1] + segments [[n_segments]]$yfr <- edge_pts$y [n_segments - 1] + + # Calculate distances and update weights for each segment + for (s in seq_len (n_segments)) { + # Calculate distance for this segment using geodesic distance + segment_xy <- data.frame ( + x = c (segments [[s]]$xfr, segments [[s]]$xto), + y = c (segments [[s]]$yfr, segments [[s]]$yto) + ) + segment_dist <- geodist::geodist (segment_xy, measure = "geodesic") [1, 2] + + # Update segment properties - preserve weight ratios exactly + segments [[s]]$d <- segment_dist + segments [[s]]$d_weighted <- segment_dist * orig_ratio + segments [[s]]$time <- segment_dist * (current_edge_1$time / current_edge_1$d) + segments [[s]]$time_weighted <- segment_dist * (current_edge_1$time / current_edge_1$d) * orig_time_ratio + + # Update edge_id to make it unique + segments [[s]]$edge_id <- paste0 (segments [[s]]$edge_id, "_", LETTERS [s]) + segments [[s]]$highway <- unique (graph_to_add [graph_to_add$tmp_graph_index == edge_pts$index [1], "highway"]) + } + + # Combine all segments + all_segments <- rbind (all_segments, do.call (rbind, segments)) + + # Add to all edges split + all_edges_split [[length (all_edges_split) + 1L]] <- all_segments + } + + # Combine all split edges + if (length (all_edges_split) == 0) { + return (graph) # No edges were split + } + + edges_split <- do.call (rbind, all_edges_split) + + # Then match edges_split back on to original graph: + graph_to_add <- graph_to_add [edges_split$n, ] + gr_cols <- gr_cols [which (!is.na (gr_cols))] + for (g in seq_along (gr_cols)) { + graph_to_add [, gr_cols [g]] <- edges_split [[names (gr_cols) [g]]] + } + graph_to_add [, "highway"] <- edges_split [["highway"]] + + # Combine original graph with new edges + result_graph <- rbind (graph, graph_to_add) + + # Update component IDs if requested + if (replace_component) { + result_graph$component <- NULL + result_graph <- dodgr::dodgr_components (result_graph) + } + result_graph$tmp_graph_index <- NULL + return (result_graph) +} diff --git a/tests/testthat/test-add-nodes-by-edge.R b/tests/testthat/test-add-nodes-by-edge.R new file mode 100644 index 00000000..c387cdb2 --- /dev/null +++ b/tests/testthat/test-add-nodes-by-edge.R @@ -0,0 +1,378 @@ +library(dplyr) +library(glue) +library(testthat) + +# test_all <- (identical (Sys.getenv ("MPADGE_LOCAL"), "true") || +# identical (Sys.getenv ("GITHUB_JOB"), "test-coverage")) +# +# skip_if (!test_all) + +library(dplyr) +dodgr_cache_off () +clear_dodgr_cache () + +test_that ("add_nodes_to_graph_by_edge single point per edge", { + + # Load a sample graph + graph <- weight_streetnet (hampi, wt_profile = "foot")%>% + mutate(graph_index=1:n()) + verts <- dodgr_vertices (graph) + + # Create a small set of points that will likely match to different edges + set.seed (1) + npts <- 5 + xy <- data.frame ( + x = min (verts$x) + runif (npts) * diff (range (verts$x)), + y = min (verts$y) + runif (npts) * diff (range (verts$y)) + ) + # Match points to graph to verify they match to different edges + pts <- match_pts_to_graph (graph, xy, distances = TRUE)%>% + mutate(xy_index=1:n())%>% + group_by(index)%>% + dplyr::filter(n()==1) + + # Filter to keep only points that match to unique edges + xy_single <- xy [pts$xy_index, ] + xy_single <- xy_single[1,] + # find bidirectional edges + bi <- graph[pts$index,]%>% + select(from_id=to_id, to_id=from_id, graph_index)%>% + inner_join(graph%>%rename(graph_index_bi=graph_index)) + bi_index <- sort(c(bi$graph_index, bi$graph_index_bi)) + + # Process with both functions + graph1 <- add_nodes_to_graph (graph, xy_single, intersections_only = TRUE, dist_tol = 0) + graph2 <- add_nodes_to_graph_by_edge (graph, xy_single, intersections_only = TRUE, dist_tol = 0) + + # Compare results + expect_equal (nrow (graph1), nrow (graph2)) + + # Compare total distance in the graph + total_dist1 <- sum (graph1$d) + total_dist2 <- sum (graph2$d) + expect_equal (total_dist1, total_dist2, tolerance = 1e-6) + + # Compare structure + expect_equal (ncol (graph1), ncol (graph2)) + + # Compare number of unique vertices + verts1 <- unique (c (graph1$from, graph1$to)) + verts2 <- unique (c (graph2$from, graph2$to)) + expect_equal (length (verts1), length (verts2)) + + + # Compare with edge to point creation + + # Process with both functions + graph1 <- add_nodes_to_graph (graph, xy_single, intersections_only = FALSE, dist_tol = 0) + graph2 <- add_nodes_to_graph_by_edge (graph, xy_single, intersections_only = FALSE, dist_tol = 0) + expect_equal (sum(graph1$d), sum(graph2$d), tolerance = 1e-6) + expect_equal (sum(graph1$time), sum(graph2$time), tolerance = 1e-6) + expect_equal (sum(graph1$time_weighted), sum(graph2$time_weighted), tolerance = 1e-6) + expect_equal (sum(graph1$d_weighted), sum(graph2$d_weighted), tolerance = 1e-6) + + + + graph2 <- add_nodes_to_graph_by_edge (graph, xy_single, intersections_only = FALSE, dist_tol = 0, wt_profile = "foot", highway="residential") + expect_equal (sum(graph1$d), sum(graph2$d), tolerance = 1e-6) + expect_equal (sum(graph1$time), sum(graph2$time), tolerance = 1e-6) + expect_equal (sum(graph1$time_weighted), sum(graph2$time_weighted), tolerance = 1e-6) + expect_equal (sum(graph1$d_weighted), sum(graph2$d_weighted), tolerance = 1e-6) + + + + graph2 <- add_nodes_to_graph_by_edge (graph, xy_single, intersections_only = FALSE, dist_tol = 0, wt_profile = "motorcar", highway="residential") + expect_equal (sum(graph1$d), sum(graph2$d), tolerance = 1e-6) + expect_gt (sum(graph1$time), sum(graph2$time)) + expect_gt (sum(graph1$time_weighted), sum(graph2$time_weighted)) +}) + +test_that ("add_nodes_to_graph_by_edge multiple points per edge", { + + # Load a sample graph + graph <- weight_streetnet (hampi, wt_profile = "foot") + verts <- dodgr_vertices (graph) + + # Create a set of points where multiple points will match to the same edge + set.seed (2) + npts <- 20 # More points increases chance of multiple matches + xy <- data.frame ( + x = min (verts$x) + runif (npts) * diff (range (verts$x)), + y = min (verts$y) + runif (npts) * diff (range (verts$y)) + ) + + # Match points to graph + pts <- match_pts_to_graph (graph, xy, distances = TRUE) + edge_counts <- table (pts$index) + + # Identify edges with multiple points + multi_point_edges <- as.integer (names (edge_counts [edge_counts > 1])) + + # Verify we have at least one edge with multiple points + expect_true (length (multi_point_edges) > 0) + + # Create a dataset with multiple points per edge + multi_point_indices <- which (pts$index %in% multi_point_edges) + xy_multi <- xy [multi_point_indices, ] + # Process with both functions + graph1 <- add_nodes_to_graph (graph, xy_multi) + graph2 <- add_nodes_to_graph_by_edge (graph, xy_multi) + + # Compare results + + # The edge-based approach should be more efficient with multiple points + # So the total distance should be less or equal + total_dist1 <- sum (graph1$d) + total_dist2 <- sum (graph2$d) + + # The edge-based approach should create fewer edges + expect_true (nrow (graph2) <= nrow (graph1)) + + # The total distance should be less or equal for the edge-based approach + expect_true (total_dist2 <= total_dist1 ) + + # Print the difference for diagnostic purposes + cat ("\nMultiple points per edge comparison:\n") + cat ("Original graph edges:", nrow (graph), "\n") + cat ("add_nodes_to_graph edges:", nrow (graph1), "\n") + cat ("add_nodes_to_graph_by_edge edges:", nrow (graph2), "\n") + cat ("add_nodes_to_graph total distance:", total_dist1, "\n") + cat ("add_nodes_to_graph_by_edge total distance:", total_dist2, "\n") + cat ("Distance ratio (edge/point):", total_dist2 / total_dist1, "\n") +}) + +test_that ("add_nodes_to_graph_by_edge with mixed point distribution", { + + # Load a sample graph + graph <- weight_streetnet (hampi, wt_profile = "foot")%>% + mutate(edge_id=as.character(edge_id))%>% + std_graph() + verts <- dodgr_vertices (graph) + + # Create a larger set of points with mixed distribution + set.seed (3) + npts <- 3 + xy <- data.frame ( + x = min (verts$x) + runif (npts) * diff (range (verts$x)), + y = min (verts$y) + runif (npts) * diff (range (verts$y)) + ) + #xy <- rbind(xy,xy) + # Process with both functions + graph1 <- add_nodes_to_graph (graph, xy, dist_tol = 0) + graph2 <- add_nodes_to_graph_by_edge (graph, xy, dist_tol = 0) + + g1 <- graph1%>% + anti_join(graph)%>% + filter(grepl("_", edge_id))%>% + tidyr::separate(edge_id, c("edge_id", "edge_id_seq"), sep="_")%>% + left_join(graph%>%select(edge_id, highway), by="edge_id")%>% + filter(highway.x!=highway.y) + + g1 <- graph1%>%anti_join(graph)%>%std_graph() + g2 <- graph2%>%anti_join(graph)%>%std_graph() + g1%>%anti_join(g2, by=c("from_lat", "from_lon", "to_lat", "to_lon", "d", "d_weighted", "highway", "time", "time_weighted"))%>%left_join(g2, by=c("from_lat", "from_lon", "to_lat", "to_lon"))%>%select(sort(names(.)))%>%View() + + # Compare results + + # The edge-based approach should generally be more efficient + total_dist1 <- sum (graph1$d) + total_dist2 <- sum (graph2$d) + cat ("Distance ratio (edge/point):", total_dist2 / total_dist1, "\n") + + # Print the difference for diagnostic purposes + cat ("\nMixed point distribution comparison:\n") + cat ("Original graph edges:", nrow (graph), "\n") + cat ("add_nodes_to_graph edges:", nrow (graph1), "\n") + cat ("add_nodes_to_graph_by_edge edges:", nrow (graph2), "\n") + cat ("add_nodes_to_graph total distance:", total_dist1, "\n") + cat ("add_nodes_to_graph_by_edge total distance:", total_dist2, "\n") + cat ("Distance ratio (edge/point):", total_dist2 / total_dist1, "\n") + + # Calculate efficiency metrics + edge_increase1 <- nrow (graph1) - nrow (graph) + edge_increase2 <- nrow (graph2) - nrow (graph) + + cat ("Edge increase (point method):", edge_increase1, "\n") + cat ("Edge increase (edge method):", edge_increase2, "\n") + cat ("Edge efficiency ratio:", edge_increase2 / edge_increase1, "\n") + + # The edge-based approach should be more efficient in terms of edges added + expect_true (edge_increase2 <= edge_increase1) + + # Verify that both graphs have the same number of unique vertices + # (excluding the intermediate vertices created during edge splitting) + verts1 <- unique (c (graph1$from_id, graph1$to_id)) + verts2 <- unique (c (graph2$from_id, graph2$to_id)) + + # The number of vertices might differ slightly due to different splitting approaches + # but the difference should be small relative to the total + vertex_diff_ratio <- abs (length (verts1) - length (verts2)) / length (verts1) + expect_true (vertex_diff_ratio < 0.1) # Allow up to 10% difference +}) + +test_that ("add_nodes_to_graph_by_edge preserves edge properties", { + + # Load a sample graph + graph <- weight_streetnet (hampi, wt_profile = "foot") + + # Create a small set of test points + set.seed (4) + npts <- 10 + verts <- dodgr_vertices (graph) + xy <- data.frame ( + x = min (verts$x) + runif (npts) * diff (range (verts$x)), + y = min (verts$y) + runif (npts) * diff (range (verts$y)) + ) + + # First, identify which edges will be split by finding the matching edges for each point + pts <- match_pts_to_graph (graph, xy, distances = TRUE) + + # Get the original edges that will be split + edges_to_split <- graph[pts$index, ] + + # Calculate the ratios for these specific edges + orig_d_ratios <- edges_to_split$d_weighted / edges_to_split$d + orig_time_ratios <- edges_to_split$time_weighted / edges_to_split$time + + # Process with both functions + graph1 <- add_nodes_to_graph (graph, xy) + graph2 <- add_nodes_to_graph_by_edge (graph, xy) + + # Check that all required columns are preserved + expect_equal (sort (names (graph)), sort (names (graph1))) + expect_equal (sort (names (graph)), sort (names (graph2))) + + # For each original edge that was split, find the corresponding new edges + for (i in seq_len(nrow(edges_to_split))) { + edge_id <- edges_to_split$edge_id[i] + + # Find new edges in graph1 that replaced this edge + # These will have edge_ids that start with the original edge_id followed by "_" + new_edges1 <- graph1[grep(paste0("^", edge_id, "_"), graph1$edge_id), ] + + # Find new edges in graph2 that replaced this edge + new_edges2 <- graph2[grep(paste0("^", edge_id, "_"), graph2$edge_id), ] + + # Skip if no matching edges found (could happen if the point was very close to a vertex) + if (nrow(new_edges1) == 0 || nrow(new_edges2) == 0) next + + # Calculate the ratios for the new edges + new_d_ratios1 <- new_edges1$d_weighted / new_edges1$d + new_time_ratios1 <- new_edges1$time_weighted / new_edges1$time + + new_d_ratios2 <- new_edges2$d_weighted / new_edges2$d + new_time_ratios2 <- new_edges2$time_weighted / new_edges2$time + + # The ratios for the new edges should be similar to the original edge + # Use mean to account for small variations due to floating point arithmetic + expect_equal (orig_d_ratios[i], mean(new_d_ratios1), tolerance = 0.01) + expect_equal (orig_time_ratios[i], mean(new_time_ratios1), tolerance = 0.01) + + expect_equal (orig_d_ratios[i], mean(new_d_ratios2), tolerance = 0.01) + expect_equal (orig_time_ratios[i], mean(new_time_ratios2), tolerance = 0.01) + + # Also check that both functions produce similar results + expect_equal (mean(new_d_ratios1), mean(new_d_ratios2), tolerance = 0.01) + expect_equal (mean(new_time_ratios1), mean(new_time_ratios2), tolerance = 0.01) + } + + # Print some diagnostic information + cat("\nEdge property preservation test:\n") + cat("Original edges to split:", nrow(edges_to_split), "\n") + cat("Original d_weighted/d ratios:", mean(orig_d_ratios), "\n") + cat("Original time_weighted/time ratios:", mean(orig_time_ratios), "\n") +}) + +test_that ("add_nodes_to_graph_by_edge handles dist_tol parameter correctly", { + + # Load a sample graph + graph <- weight_streetnet (hampi, wt_profile = "foot") + verts <- dodgr_vertices (graph) + + # Create a set of points that will be placed very close to vertices + # to test the dist_tol parameter + set.seed (5) + npts <- 10 + + # Get some random vertices from the graph + sample_verts <- verts[sample(nrow(verts), npts), ] + + # Create points that are very close to these vertices (within 1e-5 units) + xy <- data.frame( + x = sample_verts$x + rnorm(npts, 0, 1e-5), + y = sample_verts$y + rnorm(npts, 0, 1e-5) + ) + + # Test with different dist_tol values + + # With a small tolerance, points should be treated as separate + small_tol <- 1e-6 + graph1_small <- add_nodes_to_graph(graph, xy, dist_tol = small_tol) + graph2_small <- add_nodes_to_graph_by_edge(graph, xy, dist_tol = small_tol) + + # With a larger tolerance, points should be merged with existing vertices + large_tol <- 1 + graph1_large <- add_nodes_to_graph(graph, xy, dist_tol = large_tol) + graph2_large <- add_nodes_to_graph_by_edge(graph, xy, dist_tol = large_tol) + + # Compare results + + # With small tolerance, both functions should add more edges + expect_true(nrow(graph1_small) > nrow(graph)) + expect_true(nrow(graph2_small) > nrow(graph)) + + # With large tolerance, both functions should add fewer edges + # compared to the small tolerance case + expect_true(nrow(graph1_large) <= nrow(graph1_small)) + expect_true(nrow(graph2_large) <= nrow(graph2_small)) + + # The edge-based approach should be consistent with the original function + # in how it handles the tolerance parameter + small_diff_ratio <- abs(nrow(graph1_small) - nrow(graph2_small)) / nrow(graph1_small) + large_diff_ratio <- abs(nrow(graph1_large) - nrow(graph2_large)) / nrow(graph1_large) + + # The difference between the two approaches should be similar regardless of tolerance + expect_true(small_diff_ratio < 0.2) # Allow up to 20% difference + expect_true(large_diff_ratio < 0.2) + + # Print diagnostic information + cat("\ndist_tol parameter comparison:\n") + cat("Original graph edges:", nrow(graph), "\n") + cat(glue::glue("Small tolerance {small_tol}:\n")) + cat(" add_nodes_to_graph edges:", nrow(graph1_small), "\n") + cat(" add_nodes_to_graph_by_edge edges:", nrow(graph2_small), "\n") + cat(glue::glue("Large tolerance {large_tol}:\n")) + cat(" add_nodes_to_graph edges:", nrow(graph1_large), "\n") + cat(" add_nodes_to_graph_by_edge edges:", nrow(graph2_large), "\n") + cat("Edge difference ratio (small tolerance):", small_diff_ratio, "\n") + cat("Edge difference ratio (large tolerance):", large_diff_ratio, "\n") +}) + +test_that("add_nodes_to_graph_by_edge maintains consistent coordinates by ID", { + + # Load a sample graph + graph <- weight_streetnet(hampi, wt_profile = "foot") + verts <- dodgr_vertices(graph) + + # Create a set of points + set.seed(42) + npts <- 2 + xy <- data.frame( + x = min(verts$x) + runif(npts) * diff(range(verts$x)), + y = min(verts$y) + runif(npts) * diff(range(verts$y)) + ) + + # Create custom IDs for the points + xy_id <- paste0("point_", seq_len(nrow(xy))) + + # Process with add_nodes_to_graph_by_edge using custom IDs + graph_with_nodes <- add_nodes_to_graph_by_edge(graph, xy, xy_id = xy_id) + + # Get coordinates for this ID + inconsistent_coords <- graph_with_nodes %>% + dplyr::reframe(id=c(from_id, to_id), lon=c(from_lon, to_lon), lat=c(from_lat, to_lat))%>% + distinct(id, lon, lat)%>% + count(id)%>% + filter(n>1) + expect_equal(nrow(inconsistent_coords),0) +}) From 990b69f79c4a5e2fe5ffdbbc0df3413ec4f5e8d6 Mon Sep 17 00:00:00 2001 From: Eduardo Leoni Date: Fri, 4 Apr 2025 11:29:01 -0400 Subject: [PATCH 2/4] use "measure" to calculate distance, fix test, add help page --- NAMESPACE | 1 + R/match-points-by-edge.R | 10 ++-- man/add_nodes_to_graph_by_edge.Rd | 61 +++++++++++++++++++++++++ tests/testthat/test-add-nodes-by-edge.R | 38 ++++++++++----- 4 files changed, 95 insertions(+), 15 deletions(-) create mode 100644 man/add_nodes_to_graph_by_edge.Rd diff --git a/NAMESPACE b/NAMESPACE index 5b5707d9..fa3760d7 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -7,6 +7,7 @@ S3method(weight_streetnet,sc) S3method(weight_streetnet,sf) export("%>%") export(add_nodes_to_graph) +export(add_nodes_to_graph_by_edge) export(clear_dodgr_cache) export(compare_heaps) export(dodgr_cache_off) diff --git a/R/match-points-by-edge.R b/R/match-points-by-edge.R index 86b5a02b..1215905b 100644 --- a/R/match-points-by-edge.R +++ b/R/match-points-by-edge.R @@ -35,6 +35,8 @@ add_nodes_to_graph_by_edge <- function (graph, max_distance = Inf, replace_component = TRUE) { + measure <- get_geodist_measure (graph) + genhash <- function (len = 10) { paste0 (sample (c (0:9, letters, LETTERS), size = len), collapse = "") } @@ -177,13 +179,13 @@ add_nodes_to_graph_by_edge <- function (graph, new_edges$xto [1] <- new_edges$xfr [2] <- edge_pts$x [p] new_edges$yto [1] <- new_edges$yfr [2] <- edge_pts$y [p] - # Calculate distance using geodesic distance + # Calculate distance d_i <- geodist::geodist ( data.frame ( x = c (new_edges$xfr [1], new_edges$xto [1]), y = c (new_edges$yfr [1], new_edges$yto [1]) ), - measure = "geodesic" + measure = measure ) [1, 2] # Skip if distance smaller than dist_tol @@ -272,12 +274,12 @@ add_nodes_to_graph_by_edge <- function (graph, # Calculate distances and update weights for each segment for (s in seq_len (n_segments)) { - # Calculate distance for this segment using geodesic distance + # Calculate distance for this segment segment_xy <- data.frame ( x = c (segments [[s]]$xfr, segments [[s]]$xto), y = c (segments [[s]]$yfr, segments [[s]]$yto) ) - segment_dist <- geodist::geodist (segment_xy, measure = "geodesic") [1, 2] + segment_dist <- geodist::geodist (segment_xy, measure = measure) [1, 2] # Update segment properties - preserve weight ratios exactly segments [[s]]$d <- segment_dist diff --git a/man/add_nodes_to_graph_by_edge.Rd b/man/add_nodes_to_graph_by_edge.Rd new file mode 100644 index 00000000..64a4762d --- /dev/null +++ b/man/add_nodes_to_graph_by_edge.Rd @@ -0,0 +1,61 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/match-points-by-edge.R +\name{add_nodes_to_graph_by_edge} +\alias{add_nodes_to_graph_by_edge} +\title{Add nodes to a graph by matching them to the nearest edges, processing +multiple points on the same edge efficiently.} +\usage{ +add_nodes_to_graph_by_edge( + graph, + xy, + dist_tol = 1e-06, + intersections_only = FALSE, + xy_id = NULL, + wt_profile = NULL, + wt_profile_file = NULL, + highway = NULL, + max_distance = Inf, + replace_component = TRUE +) +} +\arguments{ +\item{graph}{A \code{dodgr} graph with spatial coordinates, such as a +\code{dodgr_streetnet} object.} + +\item{xy}{coordinates of points to be matched to the vertices, either as +matrix or \pkg{sf}-formatted \code{data.frame}.} + +\item{dist_tol}{Only insert new nodes if they are further from existing nodes +than this distance, expressed in units of the distance column of \code{graph}.} + +\item{intersections_only}{If \code{FALSE}} + +\item{xy_id}{Optional vector of IDs to assign to the new nodes. If NULL, IDs +will be generated automatically.} + +\item{wt_profile}{Character string specifying the weight profile to use for +edges connecting to original points (e.g., "foot", "bike", "car").} + +\item{wt_profile_file}{Character string specifying the file path to a custom +weight profile CSV file. Used only if \code{wt_profile} is NULL.} + +\item{highway}{Character string specifying the highway type for edges connecting +to original points. Must be a type defined in the weight profile.} + +\item{max_distance}{Numeric value specifying the maximum allowed distance for edges +connecting to original points. Points beyond this distance won't be connected.} + +\item{replace_component}{Logical value specifying whether to recalculate graph +component IDs after adding the new nodes and edges.} +} +\value{ +Modified version of graph with nodes added at specified locations, +original edges split at intersection points, and (unless \code{intersections_only=TRUE}) +additional edges connecting to the original points. +} +\description{ +This function extends \link{add_nodes_to_graph} by handling cases where multiple +points match to the same edge more efficiently. It splits edges at points of +perpendicular intersection and optionally creates edges connecting to the original +points. +} diff --git a/tests/testthat/test-add-nodes-by-edge.R b/tests/testthat/test-add-nodes-by-edge.R index c387cdb2..8d327818 100644 --- a/tests/testthat/test-add-nodes-by-edge.R +++ b/tests/testthat/test-add-nodes-by-edge.R @@ -20,7 +20,7 @@ test_that ("add_nodes_to_graph_by_edge single point per edge", { # Create a small set of points that will likely match to different edges set.seed (1) - npts <- 5 + npts <- 50 xy <- data.frame ( x = min (verts$x) + runif (npts) * diff (range (verts$x)), y = min (verts$y) + runif (npts) * diff (range (verts$y)) @@ -48,9 +48,11 @@ test_that ("add_nodes_to_graph_by_edge single point per edge", { expect_equal (nrow (graph1), nrow (graph2)) # Compare total distance in the graph + total_dist0 <- sum (graph$d) total_dist1 <- sum (graph1$d) total_dist2 <- sum (graph2$d) expect_equal (total_dist1, total_dist2, tolerance = 1e-6) + expect_equal (total_dist2, total_dist0, tolerance = 1e-6) # Compare structure expect_equal (ncol (graph1), ncol (graph2)) @@ -71,7 +73,9 @@ test_that ("add_nodes_to_graph_by_edge single point per edge", { expect_equal (sum(graph1$time_weighted), sum(graph2$time_weighted), tolerance = 1e-6) expect_equal (sum(graph1$d_weighted), sum(graph2$d_weighted), tolerance = 1e-6) + graph <- graph%>%filter(highway=="residential") + graph1 <- add_nodes_to_graph_by_edge (graph, xy_single, intersections_only = FALSE, dist_tol = 0) graph2 <- add_nodes_to_graph_by_edge (graph, xy_single, intersections_only = FALSE, dist_tol = 0, wt_profile = "foot", highway="residential") expect_equal (sum(graph1$d), sum(graph2$d), tolerance = 1e-6) @@ -145,32 +149,44 @@ test_that ("add_nodes_to_graph_by_edge with mixed point distribution", { # Load a sample graph graph <- weight_streetnet (hampi, wt_profile = "foot")%>% - mutate(edge_id=as.character(edge_id))%>% - std_graph() + mutate(edge_id=as.character(edge_id)) verts <- dodgr_vertices (graph) # Create a larger set of points with mixed distribution set.seed (3) - npts <- 3 + npts <- 100 xy <- data.frame ( x = min (verts$x) + runif (npts) * diff (range (verts$x)), y = min (verts$y) + runif (npts) * diff (range (verts$y)) ) #xy <- rbind(xy,xy) # Process with both functions - graph1 <- add_nodes_to_graph (graph, xy, dist_tol = 0) - graph2 <- add_nodes_to_graph_by_edge (graph, xy, dist_tol = 0) + graph1 <- add_nodes_to_graph (graph, xy, dist_tol = 0, intersections_only = TRUE) + graph2 <- add_nodes_to_graph_by_edge (graph, xy, dist_tol = 0, intersections_only = TRUE) g1 <- graph1%>% anti_join(graph)%>% filter(grepl("_", edge_id))%>% tidyr::separate(edge_id, c("edge_id", "edge_id_seq"), sep="_")%>% - left_join(graph%>%select(edge_id, highway), by="edge_id")%>% - filter(highway.x!=highway.y) + group_by(edge_id)%>% + summarise(across(c("d", "d_weighted", "time", "time_weighted"), sum))%>% + left_join(graph, by=c("edge_id"))%>% + select(sort(names(.))) - g1 <- graph1%>%anti_join(graph)%>%std_graph() - g2 <- graph2%>%anti_join(graph)%>%std_graph() - g1%>%anti_join(g2, by=c("from_lat", "from_lon", "to_lat", "to_lon", "d", "d_weighted", "highway", "time", "time_weighted"))%>%left_join(g2, by=c("from_lat", "from_lon", "to_lat", "to_lon"))%>%select(sort(names(.)))%>%View() + g2 <- graph2%>% + anti_join(graph)%>% + filter(grepl("^[0-9]+_", edge_id))%>% + tidyr::separate(edge_id, c("edge_id", "edge_id_seq"), sep="_")%>% + group_by(edge_id)%>% + summarise(across(c("d", "d_weighted", "time", "time_weighted"), sum))%>% + left_join(graph, by=c("edge_id"))%>% + select(sort(names(.))) + + #g1 <- graph1%>%anti_join(graph) + #g2 <- graph2%>%anti_join(graph) + #g1%>%anti_join(g2, by=c("from_lat", "from_lon", "to_lat", "to_lon", "d", "d_weighted", "highway", "time", "time_weighted"))%>%left_join(g2, by=c("from_lat", "from_lon", "to_lat", "to_lon"))%>%select(sort(names(.)))%>%View() + + #g2%>%anti_join(g1, by=c("from_lat", "from_lon", "to_lat", "to_lon", "d", "d_weighted", "highway", "time", "time_weighted"))%>%left_join(g1, by=c("from_lat", "from_lon", "to_lat", "to_lon"))%>%select(sort(names(.)))%>%View() # Compare results From d19ea03bb3b3b508abbf16a52274764fabb5e39f Mon Sep 17 00:00:00 2001 From: Eduardo Leoni Date: Fri, 4 Apr 2025 15:04:04 -0400 Subject: [PATCH 3/4] remove dodgr:: and dodgr::: from inside functions --- R/match-points-by-edge.R | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/R/match-points-by-edge.R b/R/match-points-by-edge.R index 1215905b..0f13ce7f 100644 --- a/R/match-points-by-edge.R +++ b/R/match-points-by-edge.R @@ -201,7 +201,7 @@ add_nodes_to_graph_by_edge <- function (graph, # Apply custom weight profile if provided if (!is.null (wt_profile) || !is.null (wt_profile_file)) { # Get weight profile - wp <- dodgr:::get_profile (wt_profile = wt_profile, file = wt_profile_file) + wp <- get_profile (wt_profile = wt_profile, file = wt_profile_file) way_wt <- wp$value [wp$way == highway] if (length (way_wt) == 0) { @@ -213,9 +213,9 @@ add_nodes_to_graph_by_edge <- function (graph, new_edges$d_weighted <- d_i / way_wt new_edges$highway <- highway # Apply additional weighting functions - new_edges <- dodgr:::set_maxspeed (new_edges, wt_profile, wt_profile_file) |> - dodgr:::weight_by_num_lanes (wt_profile) |> - dodgr:::calc_edge_time (wt_profile) + new_edges <- set_maxspeed (new_edges, wt_profile, wt_profile_file) |> + weight_by_num_lanes (wt_profile) |> + calc_edge_time (wt_profile) } else { # Use original edge's weight ratios @@ -320,7 +320,7 @@ add_nodes_to_graph_by_edge <- function (graph, # Update component IDs if requested if (replace_component) { result_graph$component <- NULL - result_graph <- dodgr::dodgr_components (result_graph) + result_graph <- dodgr_components (result_graph) } result_graph$tmp_graph_index <- NULL return (result_graph) From 2a96cda501615f5baa24982f1d1e190575b41d77 Mon Sep 17 00:00:00 2001 From: Eduardo Leoni Date: Sun, 6 Apr 2025 12:24:17 -0400 Subject: [PATCH 4/4] 1) return graph unmodified if not points left to process after cheching max distance 2) fix to ids in middle segments --- R/match-points-by-edge.R | 42 +++++++++++++++++----------------------- 1 file changed, 18 insertions(+), 24 deletions(-) diff --git a/R/match-points-by-edge.R b/R/match-points-by-edge.R index 0f13ce7f..79b67a47 100644 --- a/R/match-points-by-edge.R +++ b/R/match-points-by-edge.R @@ -35,8 +35,6 @@ add_nodes_to_graph_by_edge <- function (graph, max_distance = Inf, replace_component = TRUE) { - measure <- get_geodist_measure (graph) - genhash <- function (len = 10) { paste0 (sample (c (0:9, letters, LETTERS), size = len), collapse = "") } @@ -61,10 +59,6 @@ add_nodes_to_graph_by_edge <- function (graph, pts$xy_id <- xy_id } - # Return original graph if no points match - if (nrow (pts) == 0) { - return (graph) - } graph$tmp_graph_index <- 1:nrow (graph) @@ -93,9 +87,14 @@ add_nodes_to_graph_by_edge <- function (graph, by.y = c ("to", "from"), # Swap the order in the merge operation by.x = c ("from", "to") ) - + # Combine with original points pts <- unique (rbind (pts, pts_bi [names (pts)])) + pts <- pts[abs(pts$d_signed)<=max_distance,] # Skip pts if distance exceeds maximum allowed + # Return original graph if no points match + if (nrow (pts) == 0) { + return (graph) + } # Extract edges that need to be split edges_to_split <- graph_std [pts$index, ] @@ -179,13 +178,13 @@ add_nodes_to_graph_by_edge <- function (graph, new_edges$xto [1] <- new_edges$xfr [2] <- edge_pts$x [p] new_edges$yto [1] <- new_edges$yfr [2] <- edge_pts$y [p] - # Calculate distance + # Calculate distance using geodesic distance d_i <- geodist::geodist ( data.frame ( x = c (new_edges$xfr [1], new_edges$xto [1]), y = c (new_edges$yfr [1], new_edges$yto [1]) ), - measure = measure + measure = "geodesic" ) [1, 2] # Skip if distance smaller than dist_tol @@ -193,15 +192,10 @@ add_nodes_to_graph_by_edge <- function (graph, next } - # Skip if distance exceeds maximum allowed - if (d_i > max_distance) { - next - } - # Apply custom weight profile if provided if (!is.null (wt_profile) || !is.null (wt_profile_file)) { # Get weight profile - wp <- get_profile (wt_profile = wt_profile, file = wt_profile_file) + wp <- dodgr:::get_profile (wt_profile = wt_profile, file = wt_profile_file) way_wt <- wp$value [wp$way == highway] if (length (way_wt) == 0) { @@ -213,9 +207,9 @@ add_nodes_to_graph_by_edge <- function (graph, new_edges$d_weighted <- d_i / way_wt new_edges$highway <- highway # Apply additional weighting functions - new_edges <- set_maxspeed (new_edges, wt_profile, wt_profile_file) |> - weight_by_num_lanes (wt_profile) |> - calc_edge_time (wt_profile) + new_edges <- dodgr:::set_maxspeed (new_edges, wt_profile, wt_profile_file) |> + dodgr:::weight_by_num_lanes (wt_profile) |> + dodgr:::calc_edge_time (wt_profile) } else { # Use original edge's weight ratios @@ -254,7 +248,7 @@ add_nodes_to_graph_by_edge <- function (graph, for (s in seq_len (n_segments - 2)) { segments [[s + 1]] <- current_edge_1 segments [[s + 1]]$from <- segments [[s]]$to - segments [[s + 1]]$to <- edge_pts$proj_id [s] + segments [[s + 1]]$to <- edge_pts$proj_id [s+1] segments [[s + 1]]$xfr <- edge_pts$x [s] segments [[s + 1]]$yfr <- edge_pts$y [s] segments [[s + 1]]$xto <- edge_pts$x [s + 1] @@ -274,12 +268,12 @@ add_nodes_to_graph_by_edge <- function (graph, # Calculate distances and update weights for each segment for (s in seq_len (n_segments)) { - # Calculate distance for this segment + # Calculate distance for this segment using geodesic distance segment_xy <- data.frame ( x = c (segments [[s]]$xfr, segments [[s]]$xto), y = c (segments [[s]]$yfr, segments [[s]]$yto) ) - segment_dist <- geodist::geodist (segment_xy, measure = measure) [1, 2] + segment_dist <- geodist::geodist (segment_xy, measure = "geodesic") [1, 2] # Update segment properties - preserve weight ratios exactly segments [[s]]$d <- segment_dist @@ -289,7 +283,7 @@ add_nodes_to_graph_by_edge <- function (graph, # Update edge_id to make it unique segments [[s]]$edge_id <- paste0 (segments [[s]]$edge_id, "_", LETTERS [s]) - segments [[s]]$highway <- unique (graph_to_add [graph_to_add$tmp_graph_index == edge_pts$index [1], "highway"]) + segments [[s]]$highway <- unlist(unique (graph_to_add [graph_to_add$tmp_graph_index == edge_pts$index [1], "highway"])) } # Combine all segments @@ -305,7 +299,7 @@ add_nodes_to_graph_by_edge <- function (graph, } edges_split <- do.call (rbind, all_edges_split) - + # Then match edges_split back on to original graph: graph_to_add <- graph_to_add [edges_split$n, ] gr_cols <- gr_cols [which (!is.na (gr_cols))] @@ -320,7 +314,7 @@ add_nodes_to_graph_by_edge <- function (graph, # Update component IDs if requested if (replace_component) { result_graph$component <- NULL - result_graph <- dodgr_components (result_graph) + result_graph <- dodgr::dodgr_components (result_graph) } result_graph$tmp_graph_index <- NULL return (result_graph)