match points to graph using line intersection #103
Comments
|
Cross ref this additional |
|
Relevant Stack Overflow discussion here. |
|
Those commits add the required functionality, and introduce a breaking change: Only remaining task is to add another function to insert new points/vertices in the network, as with ropensci/stplanr/issues/237 |
library (dodgr)
packageVersion ("dodgr")
#> [1] '0.2.14.90'
graph <- weight_streetnet (hampi, wt_profile = "foot")
dim (graph)
#> [1] 6813 15
verts <- dodgr_vertices (graph)
x11 ()
par (mfrow = c (1, 2))
plot (NULL, xlim = range (verts$x), ylim = range (verts$y))
segments (graph$from_lon, graph$from_lat, graph$to_lon, graph$to_lat)
set.seed (2)
npts <- 10
xy <- data.frame (
x = min (verts$x) + runif (npts) * diff (range (verts$x)),
y = min (verts$y) + runif (npts) * diff (range (verts$y))
)
points (xy$x, xy$y, pch = 19, col = "red")
graph <- add_nodes_to_graph (graph, xy)
dim (graph)
#> [1] 6833 15
plot (NULL, xlim = range (verts$x), ylim = range (verts$y))
segments (graph$from_lon, graph$from_lat, graph$to_lon, graph$to_lat)
points (xy$x, xy$y, pch = 19, col = "red")Created on 2022-08-11 by the reprex package (v2.0.1) Adding 10 points creates 20 new edges, because each one is duplicated in bi-directional form. Concluding notesRoutines currently presume street networks with EPSG4326 only; the matching-by-perpendicular-intersection routines are applied directly to projected coordinates, which merely presumes non-planarity to be negligible (so presumes graphs cover relatively small areas), and final distances are geodesic. Can revisit and modify later if this routine starts being used more widely. |
|
Re-opening to add ability to return actual distances to nearest edge. These distances can be very useful for routing tasks, to construct some measure of weighted distance to some kinds of points of attraction (or whatever). |
|
Those commits add the distance to the output of the C++function. It's not yet used, but can be added directly to the R output from there. ###TODO
Not sure the best way to do that? Likely get the actual point of intersection and return the coordinates of that so that a final application in R of ˋgeodistˋ will then give the correct distance. That will require tweaking the C++ code to get the actual coordinates of the intersection point as well as the distances. |
|
That commit adds the geometric distances. Just need to run a few test cases to make sure everything is okay, and then modify the return value to include those distances. |
|
What is the current status of this issue? I was looking for something that allow me to similar function as The current The graph plotted with and without using add_nodes_to_graph() is illustrated as the first and second image below: the red points are the start and end points. The difference is that graph has two more edges (which was explained above for bi-directional) for each node I want to add. I was expecting something that looks like the third image below.
The current way to add node to a graph creates a weird effect on my routing and makes my start and end points (which are centroids of buildings) to be part of the graph (which they are not.. the graph is the road network). I realized this will work well if the vertex-to-be-added is already ON the graph (then the current mode makes perfect sense). Is there an easy way to find the closest matching point on the graph if the given point (e.g., in my case, building centroids) is not on the graph? I need to do this for many pairs of ODs. In case you need it, this is the code for making the graphics above |
|
Thanks @xiaofanliang, short answer is yes, that's all possible, but not yet sufficiently well documented. I'll get a longer answer to you sometime next week |
|
@xiaofanliang Main problem was that the function docs had not been updating properly for quite some time. Updated version of |
|
@mpadge Thanks! I will experiment with my codes. |
|
@mpadge I tried it out with an example with the latest dev branch download of Here I attached a replicable example to work with. This is what the data looks like. I would like to find the closest perpendicular points on the graph for reach red dots.
If I run |
|
@xiaofanliang I guess by "perpendicular points" you mean the actual points of intersection between the pependicular line and the closest edge? Those points of intersection are calculated internally, but not returned by the function. I guess it would be useful to have a way to access those. I could just add two extra columns when |
library (dodgr)
packageVersion ("dodgr")
#> [1] '0.2.20.22'
graph <- weight_streetnet (hampi, wt_profile = "foot")
v <- dodgr_vertices (graph)
# A random point:
set.seed (1L)
n <- 10L
xy <- apply (v [, c ("x", "y")], 2, function (i) runif (n, min = min (i), max = max (i)))
d <- match_pts_to_graph (graph, xy, distances = TRUE)
print (d)
#> index d_signed x y
#> 1 3395 0.00000 76.42341 15.31717
#> 2 3395 0.00000 76.42341 15.31717
#> 3 5271 0.00000 76.44288 15.34326
#> 4 6325 0.00000 76.48064 15.32658
#> 5 6225 -623.02542 76.39593 15.35226
#> 6 2469 0.00000 76.47957 15.33146
#> 7 31 0.00000 76.48381 15.33999
#> 8 3247 0.00000 76.45306 15.36033
#> 9 3539 0.00000 76.44780 15.32178
#> 10 6201 54.05596 76.38011 15.34666
library (sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.4, PROJ 9.2.0; sf_use_s2() is TRUE
plot (hampi$geometry)
points (xy, col = "orange")
points (d [, c ("x", "y")], col = "red", pch = 19)
segments (xy [, 1], xy [, 2], d$x, d$y, col = "red")Created on 2023-05-08 with reprex v2.0.2 (Note that intersection lines will only look perpendicular when the plot rectangle is truly square, which is almost never will be.) |
|
@mpadge Yes! That is what I am looking for. It seems like the new updates have not been integrated in the I can explain my motivation for asking for this functionality a bit more and maybe it will help you figure out what is the best way to implement this. My start and end points for routing are building centroids. I would like to contract my graph (road networks) via |
|
@xiaofanliang That's all on main branch now, so good for you to use straight away. Thanks for explaining your general use case. I face exactly that issue a lot too, but in my case it's always sufficient just to match to vertices and ensure those are passed as the |
|
Sorry, I hadn't pushed the commit. It's now up. |
|
I experimented today with In short, I use However, when I calculate distance, the distances on the uncontracted graph and contracted graph are very different, even if I use the same starting/ending points (i.e., intersection points returned by
|
|
Yep, definitely something wrong there: library (dodgr)
graph <- weight_streetnet (hampi, wt_profile = "foot")
graph <- graph[graph$component == 1, ]
v <- dodgr_vertices (graph)
set.seed (2L)
n <- 10L
xy <- apply (v [, c ("x", "y")], 2, function (i) runif (n, min = min (i), max = max (i)))
d <- match_pts_to_graph (graph, xy, distances = TRUE)
graph2 = dodgr_contract_graph(graph)
graph2 = add_nodes_to_graph(graph2, d[,c('x','y')])
i <- c (1, 7)
dodgr_dists(graph, from = xy[i [1],], to = xy[i [2],], pairwise = T)
#> [,1]
#> [1,] 616.9445
dodgr_dists(graph2, from = xy[i [1],], to = xy[i [2],], pairwise = T)
#> [,1]
#> [1,] 5262.223
# dodgr_paths(graph, from = xy[i [1],], to = xy[i [2],], pairwise = T)
# dodgr_paths(graph2, from = xy[i [1],], to = xy[i [2],], pairwise = T)
p1 <- dodgr_paths(graph, from = xy[i [1],], to = xy[i [2],], pairwise = T)
v1 <- dodgr_vertices (graph)
p1 <- v1 [match (p1 [[1]] [[1]], v1$id), ]
p2 <- dodgr_paths(graph2, from = xy[i [1],], to = xy[i [2],], pairwise = T)
v2 <- dodgr_vertices (graph2)
p2 <- v2 [match (p2 [[1]] [[1]], v2$id), ]
plot (NULL, xlim = range (rbind (p1, p2)$x), ylim = range (rbind (p1, p2)$y))
lines (p1$x, p1$y, col = "red")
points (p1$x, p1$y, pch = 19, col = "red")
points (p1$x [1], p1$y [1], pch = 19, cex = 2, col = "red")
points (p1$x [nrow (p1)], p1$y [nrow (p1)], pch = 19, cex = 2, col = "orange")
lines (p2$x, p2$y, col = "lawngreen")
points (p2$x, p2$y, pch = 19, col = "lawngreen")
points (p2$x [1], p2$y [1], pch = 19, cex = 1.5, col = "lawngreen")
points (p2$x [nrow (p2)], p2$y [nrow (p2)], pch = 19, cex = 1.5, col = "blue")Created on 2023-05-10 with reprex v2.0.2 |
|
That example is a poor one, because the start and end points both match on to the same vertex of the original graph. But the problem definitely lies in the |
|
So the problem seems to arise when multiple points submitted to library (dodgr)
graph <- weight_streetnet (hampi, wt_profile = "foot")
graph <- graph[graph$component == 1, ]
v <- dodgr_vertices (graph)
set.seed (2L)
n <- 10L
xy <- apply (v [, c ("x", "y")], 2, function (i) runif (n, min = min (i), max = max (i)))
d <- match_pts_to_graph (graph, xy, distances = TRUE)
graph2 = dodgr_contract_graph(graph)
graph2 = add_nodes_to_graph(graph2, d[,c('x','y')])
# Important note: The match_pts_to_verts call should here
# include the intersection vertices, not the actual xy coordinates
# xy_id1 <- v$id [match_pts_to_verts (v, xy)]
xy_id1 <- v$id [match_pts_to_verts (v, d [, c ("x", "y")])]
v2 <- dodgr_vertices (graph2)
# xy_id2 <- v2$id [match_pts_to_verts (v2, xy)]
xy_id2 <- v2$id [match_pts_to_verts (v2, d [, c ("x", "y")])]
i <- c (1, 8)
# These should be the same, or almost identical, and indeed are:
v [match (xy_id1 [i [1]], v$id), ]
#> id x y component n
#> 2500 7794251384 76.43668 15.32657 1 1255
v2 [match (xy_id2 [i [1]], v2$id), ]
#> id x y component n
#> 528 OS7CwRELzG 76.43668 15.32657 1 210
# And these are also nearly the same:
v [match (xy_id1 [i [2]], v$id), ]
#> id x y component n
#> 3964 1376769158 76.47739 15.31274 1 1958
v2 [match (xy_id2 [i [2]], v2$id), ]
#> id x y component n
#> 556 mO473txWcg 76.47733 15.31263 1 228
p1 <- dodgr_paths(graph, from = xy[i [1],], to = xy[i [2],], pairwise = T)
v1 <- dodgr_vertices (graph)
p1 <- v1 [match (p1 [[1]] [[1]], v1$id), ]
p2 <- dodgr_paths(graph2, from = xy[i [1],], to = xy[i [2],], pairwise = T)
v2 <- dodgr_vertices (graph2)
p2 <- v2 [match (p2 [[1]] [[1]], v2$id), ]
plot (NULL, xlim = range (rbind (p1, p2)$x), ylim = range (rbind (p1, p2)$y))
lines (p1$x, p1$y, col = "red")
points (p1$x, p1$y, pch = 19, col = "red")
text (p1$x [1], p1$y [1], labels = "start", pos = 3, col = "red")
text (p1$x [nrow (p1)], p1$y [nrow (p1)], labels = "end", pos = 3, col = "red")
lines (p2$x, p2$y, col = "blue")
points (p2$x, p2$y, pch = 19, col = "blue")
text (p2$x [1], p2$y [1], labels = "start", pos = 3, col = "blue")
text (p2$x [nrow (p2)], p2$y [nrow (p2)], labels = "end", pos = 3, col = "blue")
points (d [, c ("x", "y")], pch = 19, col ="gray", cex = 2)Created on 2023-05-10 with reprex v2.0.2 That excursion extends out to an additional point that matches on the same edge as the start point. It is then inadvently added to the graph in the |
|
The above commit brings us here, noting that the workflow is slightly different to the above. I think this is the appropriate workflow here, which is:
library (dodgr)
packageVersion ("dodgr")
#> [1] '0.2.20.25'
graph <- weight_streetnet (hampi, wt_profile = "foot")
graph <- graph[graph$component == 1, ]
v0 <- dodgr_vertices (graph)
set.seed (2L)
n <- 10L
xy <- apply (v0 [, c ("x", "y")], 2, function (i) runif (n, min = min (i), max = max (i)))
d <- match_pts_to_graph (graph, xy, distances = TRUE)
graph = add_nodes_to_graph(graph, d[,c('x','y')])
v <- dodgr_vertices (graph)
v_new <- v$id [which (!v$id %in% v0$id)]
graph2 <- dodgr_contract_graph (graph, verts = v_new)
xy_id1 <- v$id [match_pts_to_verts (v, d [, c ("x", "y")])]
v2 <- dodgr_vertices (graph2)
xy_id2 <- v2$id [match_pts_to_verts (v2, d [, c ("x", "y")])]
i <- c (1, 8)
dodgr_dists(graph, from = xy_id1 [i [1]], to = xy_id1 [i [2]], pairwise = T)
#> [,1]
#> [1,] 8162.554
dodgr_dists(graph2, from = xy_id2 [i [1]], to = xy_id2 [i [2]], pairwise = T)
#> [,1]
#> [1,] 9085.145
p1 <- dodgr_paths(graph, from = xy_id1 [i [1]], to = xy_id1 [i [2]], pairwise = T)
v1 <- dodgr_vertices (graph)
p1 <- v1 [match (p1 [[1]] [[1]], v1$id), ]
p2 <- dodgr_paths(graph2, from = xy_id2 [i [1]], to = xy_id2 [i [2]], pairwise = T)
v2 <- dodgr_vertices (graph2)
p2 <- v2 [match (p2 [[1]] [[1]], v2$id), ]
plot (NULL, xlim = range (rbind (p1, p2)$x), ylim = range (rbind (p1, p2)$y))
segments (graph2$from_lon, graph2$from_lat, graph2$to_lon, graph2$to_lat, col = "gray")
lines (p1$x, p1$y, col = "red")
points (p1$x, p1$y, pch = 19, col = "red")
text (p1$x [1], p1$y [1], labels = "start", pos = 3, col = "red")
text (p1$x [nrow (p1)], p1$y [nrow (p1)], labels = "end", pos = 3, col = "red")
lines (p2$x, p2$y, col = "blue")
points (p2$x, p2$y, pch = 19, col = "blue")
text (p2$x [1], p2$y [1], labels = "start", pos = 3, col = "blue")
text (p2$x [nrow (p2)], p2$y [nrow (p2)], labels = "end", pos = 3, col = "blue")
points (d [, c ("x", "y")], pch = 19, col ="gray", cex = 2)
text (d [, c ("x", "y")], labels = seq_len (nrow (d)), pos = 1)Created on 2023-05-10 with reprex v2.0.2 I think that mostly fixes it. The commit introduced a That all makes sense to me, but still leaves the issue that the general workflow we want here should be able to yield identical distances from both versions of the graph, so there's still some work to be done here. |
|
That commit extends the code so that new points can be directly submitted to library (dodgr)
packageVersion ("dodgr")
#> [1] '0.2.20.26'
graph <- weight_streetnet (hampi, wt_profile = "foot")
graph <- graph[graph$component == 1, ]
v0 <- dodgr_vertices (graph)
set.seed (2L)
n <- 10L
xy <- apply (v0 [, c ("x", "y")], 2, function (i) runif (n, min = min (i), max = max (i)))
graph = add_nodes_to_graph(graph, xy)
v <- dodgr_vertices (graph)
v_new <- v$id [which (!v$id %in% v0$id)]
graph2 <- dodgr_contract_graph (graph, verts = v_new)
xy_id1 <- v$id [match_pts_to_verts (v, xy)]
v2 <- dodgr_vertices (graph2)
xy_id2 <- v2$id [match_pts_to_verts (v2, xy)]
i <- c (1, 8)
dodgr_dists(graph, from = xy_id1 [i [1]], to = xy_id1 [i [2]], pairwise = T)
#> [,1]
#> [1,] 8800.232
dodgr_dists(graph2, from = xy_id2 [i [1]], to = xy_id2 [i [2]], pairwise = T)
#> [,1]
#> [1,] 8800.232
p1 <- dodgr_paths(graph, from = xy_id1 [i [1]], to = xy_id1 [i [2]], pairwise = T)
v1 <- dodgr_vertices (graph)
p1 <- v1 [match (p1 [[1]] [[1]], v1$id), ]
p2 <- dodgr_paths(graph2, from = xy_id2 [i [1]], to = xy_id2 [i [2]], pairwise = T)
v2 <- dodgr_vertices (graph2)
p2 <- v2 [match (p2 [[1]] [[1]], v2$id), ]
plot (NULL, xlim = range (rbind (p1, p2)$x), ylim = range (rbind (p1, p2)$y))
segments (graph2$from_lon, graph2$from_lat, graph2$to_lon, graph2$to_lat, col = "gray")
lines (p1$x, p1$y, col = "red")
points (p1$x, p1$y, pch = 19, col = "red")
text (p1$x [1], p1$y [1], labels = "start", pos = 3, col = "red")
text (p1$x [nrow (p1)], p1$y [nrow (p1)], labels = "end", pos = 3, col = "red")
lines (p2$x, p2$y, col = "blue")
points (p2$x, p2$y, pch = 19, col = "blue")
text (p2$x [1], p2$y [1], labels = "start", pos = 3, col = "blue")
text (p2$x [nrow (p2)], p2$y [nrow (p2)], labels = "end", pos = 3, col = "blue")
points (xy, pch = 19, col ="gray", cex = 2)
text (xy, labels = seq_len (nrow (xy)), pos = 1)Created on 2023-05-10 with reprex v2.0.2 Now just have to expose an additional parameter, something like |
|
That commit then means the the plot immediately above is generated with graph = add_nodes_to_graph(graph, xy, intersections_only = FALSE) # defaultAnd the one prior to that, where the start points don't quite match up, and the distances differ, arises from graph = add_nodes_to_graph(graph, xy, intersections_only = TRUE)The workflow described above is then slightly simplified to one step only:
With the default graph2 <- dodgr_contract_graph (graph, verts = v_new)to graph2 <- dodgr_contract_graph (graph)and the lines to identify newly inserted vertices removed. Let me know how that goes for you @xiaofanliang, and thanks for the extensive testing 👍 |
|
Thanks for the extensive updates as well! Here is an example of the complete workflow. I would like to clarify when to use I notice that in the example below, when I set |
Current
match_points_to_graph()just matches to nearest graph vertex. This could be improved by finding the closest edge to a matched point, then choosing the closest vertex of that edge. This extends from suggestions of @JimShady in stplanr issue#331. Implementing this will require a custom line intersection routine.The text was updated successfully, but these errors were encountered: