UrbanAnalyst / dodgr

Distances on Directed Graphs in R
https://urbananalyst.github.io/dodgr/
127 stars 16 forks source link

cannot calculate dists when end/start is oneway #152

Closed JBSP-code closed 3 years ago

JBSP-code commented 3 years ago

First of all, thank you for coding dodgr, it's a blast to work with it.

I'm currently calculating drive times to retail stores in Switzerland and get systematically NAs for specific stores. I suspect that the issue occurs, whenever the nearest node to the store coordinate lies on a edge that is oneway and makes the routing therefore technically impossible, i.e. the last edge when the store is a destination is a oneway street away from the store or the starting point of the routing lies next to a node which is a oneway not leading to the rest of the street network. Below is a short example in which I defined some exemplary destination coordinates (from which the first one are the true coordinates from my dataset).

library(dodgr)

bb <- osmdata::getbb("bern")
npts <- 10
xy <- apply (bb, 1, function (i) min (i) + runif (npts) * diff (i))
destinations <- data.frame(x=c(7.455, 7.4549, 7.4552, 7.4549), 
                           y=c(46.9593, 46.9591, 46.9594, 46.9594))

net <- dodgr_streetnet_sc(bb)
net_wt <- weight_streetnet(net, wt_profile = "motorcar")
net_wt <- dodgr_contract_graph(net_wt)
net_wt <- net_wt[net_wt$component == 1, ]

dists <- dodgr_times(graph = net_wt, 
                     from = destinations,
                     to = xy, 
                     shortest = FALSE)

image

I then inspected the street network using QGIS and filtering the .osm using only streets that have a weighting profile as can be determined with dplyr::filter(weighting_profiles$weighting_profiles, name == "motorcar", !is.na(max_speed)). In the image below the center point corresponds to the first and "real" coordinates. I suspect that the center point and the one to the upper right get no drive times as they are located closest to the node on the oneway edge.

image

Do you by any chance know how I can address this issue or circumvent the problem? Your help is greatly appreciated, thank you.

mpadge commented 3 years ago

Thanks @theFippo for the encouraging words. I'm not exactly sure what the issue is, but one thing that jumps out at me is that you've submitted your destinations as from points. Reversing the order and routing to those destinations works as expected. Note that the code includes extra lines using match_points_to_graph to ensure that your routing points are included in the contracted graph. In your case, this generally makes no difference, but is good practice to include anyway to obtain more precise routing results.

library (dodgr)

bb <- osmdata::getbb("bern")
npts <- 10
xy <- apply (bb, 1, function (i) min (i) + runif (npts) * diff (i))
destinations <- data.frame(x=c(7.455, 7.4549, 7.4552, 7.4549), 
                           y=c(46.9593, 46.9591, 46.9594, 46.9594))
destinations <- destinations [c (1, 3), ] # 2 and 4 work fine

net <- dodgr_streetnet_sc(bb)
net_wt <- weight_streetnet(net, wt_profile = "motorcar")
#> Loading required namespace: geodist
#> Loading required namespace: dplyr
v <- dodgr_vertices (net_wt)
pts <- match_points_to_graph (v, rbind (xy, destinations))
net_c <- dodgr_contract_graph(net_wt, verts = pts)
net_c <- net_c[net_c$component == 1, ]

dodgr_times(graph = net_c, from = xy, to = destinations, shortest = FALSE)
#>                    1         3
#> 1782513363  443.1134  443.1134
#> 567931066  1206.4369 1206.4369
#> 3359436253 1462.7364 1462.7364
#> 303844605  1588.7159 1588.7159
#> 6594944004  958.9605  958.9605
#> 1070848269  713.3837  713.3837
#> 4396806187  781.4027  781.4027
#> 4070317559 1182.8574 1182.8574
#> 290798531   731.3850  731.3850
#> 277248573   904.4553  904.4553
dodgr_times(graph = net_c, from = destinations, to = xy, shortest = FALSE)
#>   1782513363 567931066 3359436253 303844605 6594944004 1070848269 4396806187
#> 1         NA        NA         NA        NA         NA         NA         NA
#> 3         NA        NA         NA        NA         NA         NA         NA
#>   4070317559 290798531 277248573
#> 1         NA        NA        NA
#> 3         NA        NA        NA

# Then examine each edge to and from the destination points:
pts <- match_points_to_graph (v, destinations)
ids <- v$id [pts]
message ("destinations: [", paste0 (ids, collapse = ", "), "]")
#> destinations: [597213075, 5034012341]
data.frame (net_wt [match (ids, net_wt$.vx0), ]) # from destination
#>        .vx0      .vx1      edge_   .vx0_x   .vx0_y   .vx1_x  .vx1_y        d
#> 1 597213075 565552518 6z6Qg5xYwK 7.454809 46.95922 7.454769 46.9594 20.65921
#> 2      <NA>      <NA>       <NA>       NA       NA       NA      NA       NA
#>    object_  highway oneway.bicycle oneway.bicycle.1 d_weighted     time
#> 1 46711463 tertiary           <NA>             <NA>   29.51316 1.859329
#> 2     <NA>     <NA>           <NA>             <NA>         NA       NA
#>   time_weighted component
#> 1      2.656184         1
#> 2            NA        NA
data.frame (net_wt [match (ids, net_wt$.vx1), ]) # to destination
#>         .vx0       .vx1      edge_   .vx0_x   .vx0_y   .vx1_x   .vx1_y        d
#> 1  565552518  597213075 i3Ux6Ch2F9 7.454769 46.95940 7.454809 46.95922 20.65921
#> 2 7524331429 5034012341 yaoiLKHgLo 7.454845 46.95949 7.455148 46.95953 23.30962
#>     object_  highway oneway.bicycle oneway.bicycle.1 d_weighted     time
#> 1  46711463 tertiary           <NA>             <NA>   29.51316 1.859329
#> 2 515469550  service           <NA>             <NA>   58.27405 5.594309
#>   time_weighted component
#> 1      2.656184         1
#> 2     13.985772         1

Created on 2020-12-31 by the reprex package (v0.3.0)

And there is no edge at all leading from the second of those destination points, which is why that does not work. The first point has an edge leading away, but that must somehow not connect with the rest of the network, so that both "destination" points can not be routed from when submitted as origin points. You'd have to trace the network a bit further to find out what the reason might be,

That result also shows, however, that the oneway flag is not in the network at all, even though some ways are indeed flagged as oneway for motorcar routing, so there's obviously a bug there, for which I'll open a separate issue.

JBSP-code commented 3 years ago

Thanks a lot for your help and pointing me to the right direction. Indeed, I inattentively messed up from and to locations, so thanks for that 👍 . I began by routing store locations to store locations when I noticed the issue and added those random "starting" points, which ended up in truth as destinations and the "destinations" as starting points. Anyways, by following up inspecting the matched vertices, as you have pointed them out in your code, I found that they indeed are matching to the end vertex of the oneway street, i. e. there is no way from there to get to the rest of the network. (I used match_points_to_graph() after contraction, which is why both starting points now match to the same vertex.) Here again marked in QGIS, where the left yellow marked point is the one appearing in your code but not in mine and the other one appearing in both your code as well as mine below) image To address the issue I proceeded to remove all edges of oneway streets that are starting or end points of the network, effectively reducing the available set of matching locations for the routing points .(A oneway street away from the network at the start make starting points unroutable, oneway streets at the destinations into the network make destinations unreachable. My example only shows one case, but the other appeared in other cities.)

library(dodgr)
library(magrittr)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

bb <- osmdata::getbb("bern")
npts <- 10
xy <- apply (bb, 1, function (i) min (i) + runif (npts) * diff (i))

# set some destination coordinates (of which the first one is a "real" destination)
destinations <- data.frame(x=c(7.455, 7.4549, 7.4552, 7.4549), 
                           y=c(46.9593, 46.9591, 46.9594, 46.9594))

net <- dodgr_streetnet_sc(bb)

net_wt <- weight_streetnet(net, wt_profile = "motorcar")
#> Loading required namespace: geodist
net_wt <- dodgr_contract_graph(net_wt)
net_wt <- net_wt[net_wt$component == 1, ]

dodgr_times(graph = net_wt, from = destinations, to = xy, shortest = FALSE) # 1 and 3 nonroutable
#>   538174292 5826910813 1078849513 4030496541 288961183 4173969932 1903883989
#> 1        NA         NA         NA         NA        NA         NA         NA
#> 2  752.3009   997.0716   663.9354   1125.662  723.7828   1202.399   1545.119
#> 3        NA         NA         NA         NA        NA         NA         NA
#> 4  765.3846  1010.1552   677.0191   1138.746  736.8665   1215.483   1558.202
#>   5483555833 4517774829 1258785678
#> 1         NA         NA         NA
#> 2   397.1025   860.4033   1438.032
#> 3         NA         NA         NA
#> 4   408.9318   873.4870   1451.115

destinations <- destinations [c (1, 3), ] # 2 and 4 work fine

v <- dodgr_vertices (net_wt)
pts <- match_points_to_graph (v, destinations)
ids <- v$id [pts]
message ("destinations: [", paste0 (ids, collapse = ", "), "]")
#> destinations: [5034012341, 5034012341]

# both point now match to the ending vertex of the oneway street
net_wt[net_wt$.vx1 %in% c(ids),]
#> # A tibble: 1 x 16
#>   .vx0  .vx1  edge_ .vx0_x .vx0_y .vx1_x .vx1_y     d object_ highway
#>   <chr> <chr> <chr>  <dbl>  <dbl>  <dbl>  <dbl> <dbl> <chr>   <chr>  
#> 1 5034~ 5034~ a210~   7.45   47.0   7.46   47.0  30.5 515469~ service
#> # ... with 6 more variables: oneway.bicycle <chr>, `oneway:bicycle` <chr>,
#> #   d_weighted <dbl>, time <dbl>, time_weighted <dbl>, component <int>
net_wt[net_wt$.vx0 %in% c(ids),]
#> # A tibble: 0 x 16
#> # ... with 16 variables: .vx0 <chr>, .vx1 <chr>, edge_ <chr>, .vx0_x <dbl>,
#> #   .vx0_y <dbl>, .vx1_x <dbl>, .vx1_y <dbl>, d <dbl>, object_ <chr>,
#> #   highway <chr>, oneway.bicycle <chr>, `oneway:bicycle` <chr>,
#> #   d_weighted <dbl>, time <dbl>, time_weighted <dbl>, component <int>

# as oneway streets at the end/start are causing trouble, 
# use an anti join by itself on vertices that are not leading back or further, 
# i.e. the ending vertex is itself not on an edge leading back to the starting vertex or to any other vertex

# nonroutable starting edges, i.e. oneway edges leading away from the rest
dead_starts <- unique(net_wt[c(".vx0",".vx1","edge_")]) %>% 
  as_tibble() %>%
  anti_join(., ., by = c(".vx0" = ".vx1")) %>% 
  select(edge_) %>% 
  distinct()

# nonroutable ending edges, i.e. oneway edges leading into the street network but not from endpoints back into
dead_ends <- unique(net_wt[c(".vx0",".vx1","edge_")]) %>% 
  as_tibble() %>%
  anti_join(., ., by = c(".vx1" = ".vx0")) %>% 
  select(edge_) %>% 
  distinct()

# edges to be removed
c(dead_ends[["edge_"]],dead_starts[["edge_"]])
#>  [1] "a239891"    "a239837"    "a238632"    "a235559"    "a235107"   
#>  [6] "a234991"    "a233858"    "a231326"    "a229805"    "a236685"   
#> [11] "a220585"    "a220239"    "a210185"    "a193233"    "a193181"   
#> [16] "a205505"    "a128515"    "a238404"    "a163674"    "a201222"   
#> [21] "a238117"    "a238044"    "a237445"    "a235411"    "a235106"   
#> [26] "a233827"    "a225866"    "a220495"    "a197458"    "a197162"   
#> [31] "a193219"    "a206255"    "a188581"    "a128512"    "a190909"   
#> [36] "a207338"    "a193418"    "a205330"    "a171078"    "eWFb2c7jJY"
#> [41] "19k6TrkJ80"

net_wt <- net_wt[!(net_wt$edge_ %in% c(dead_ends[["edge_"]],dead_starts[["edge_"]])),]

dodgr_times(graph = net_wt, from = destinations, to = xy, shortest = FALSE)
#>   538174292 5826910813 1078849513 4030496541 288961183 4173969932 1903883989
#> 1  765.3846   1010.155   677.0191   1138.746  736.8665   1215.483   1558.202
#> 3  708.9154   1012.913   679.7765   1141.503  739.6239   1218.240   1560.960
#>   5483555833 4517774829 1258785678
#> 1   408.9318   873.4870   1451.115
#> 3   385.2970   876.2444   1453.873

Created on 2021-01-02 by the reprex package (v0.3.0)

Of course, this approach seems to be a bit extreme, as I effectively remove streets that are in principle useable. However, I ended up only removing 41 edges out of 30953 which I consider to be ok.

As for me, the issue is taken care of. I'm not sure whether this issue occurs generally, as I've not seen it anywhere here before. Again, thank you for your help!