4

i would like to change weight of parts of a route based on gps coordinates. For that i would like to get gps coordinates of an edge of a calculated route then compare them with a list of coordinates i have and if the coordinates in my list matches the coordinates of an endge of a route, i would like to change the weight of that edge. Currently i have the code that calculates the route and changes the weight of the whole route. I get the coordinates of the route but i just cant get the steps requiered to get back to the graph.. my brain just shuts down :)

library(osmar)
library(igraph)

### Get data ----
src <- osmsource_api(url = "https://api.openstreetmap.org/api/0.6/")
muc_bbox <- center_bbox(11.575278, 48.137222, 1000, 1000)
muc <- get_osm(muc_bbox, src)

### Reduce to highways: ----
hways <- subset(muc, way_ids = find(muc, way(tags(k == "highway"))))
hways <- find(hways, way(tags(k == "name")))
hways <- find_down(muc, way(hways))
hways <- subset(muc, ids = hways)

#### Plot data ----
## Plot complete data and highways on top:
plot(muc)
plot_ways(muc, col = "lightgrey")
plot_ways(hways, col = "coral", add = TRUE)

### Define route start and end nodes: ----
id<-find(muc, node(tags(v %agrep% "Sendlinger Tor")))[1]
hway_start_node <-find_nearest_node(muc, id, way(tags(k == "highway"))) 
hway_start <- subset(muc, node(hway_start_node))

id <- find(muc, node(attrs(lon > 11.58 & lat > 48.15)))[1]
hway_end_node <- find_nearest_node(muc, id, way(tags(k == "highway")))
hway_end <- subset(muc, node(hway_end_node))

## Add the route start and and nodes to the plot:
plot_nodes(hway_start, add = TRUE, col = "red", pch = 19, cex = 2)
plot_nodes(hway_end, add = TRUE, col = "red", pch = 19, cex = 2)

### Create street graph ----
gr <- as.undirected(as_igraph(hways))

### Compute shortest route: ----
# Calculate route
route <- function(start_node,end_node) {
  get.shortest.paths(gr,
                     from = as.character(start_node),
                     to = as.character(end_node), 
                     mode = "all")[[1]][[1]]}
# Plot route
plot.route <- function(r,color) {
  r.nodes.names <- as.numeric(V(gr)[r]$name)
  r.ways <- subset(hways, ids = osmar::find_up(hways, node(r.nodes.names)))
  plot_ways(r.ways, add = TRUE, col = color, lwd = 2)
}


r <- route(hway_start_node,hway_end_node)
color <- colorRampPalette(c("springgreen","royalblue"))(nways)[numway]
plot.route(r,color)

route_nodes <- as.numeric(V(gr)[r]$name)
#We construct a newosmarobject containing only elements 
#related to the nodes defining the route:

route_ids <- find_up(hways, node(route_nodes))
route_muc <- subset(hways, ids = route_ids)

#Route details.
#In order to present route details like street names,
#distances, and directions we have to work directly on the internals of the osmar objects.
#We start by extracting the route’s node IDs (which are in the correct order) 
#and the way IDs       (whichwe have to order)
#where the nodes are members:

node_ids <- route_muc$nodes$attrs$id

way_ids <- local({
  w <- match(node_ids, route_muc$ways$refs$ref)
  route_muc$ways$refs$id[w]
})

#Then we extract the names of the ways in the correct order:>
way_names <- local({
  n <- subset(route_muc$ways$tags, k == "name")
  n[match(way_ids, n$id), "v"]
})

#The next step is to extract the nodes’ coordinates,>
node_coords <- route_muc$nodes$attrs[, c("lon", "lat")]

#and to compute the distances (meters) and the bearings (degrees) 
#between successive nodes (using thepackagegeosphere):

node_dirs <- local({
      n <- nrow(node_coords)
      from <- 1:(n - 1)
      to <- 2:n
      cbind(dist = c(0, distHaversine(node_coords[from,], node_coords[to,])),
            bear = c(0, bearing(node_coords[from,], node_coords[to,])))
    })



#Finally, we pack together all the information, 
#and additionally compute the cumulative distance

route_details <- data.frame(way_names, node_dirs)
route_details$cdist <- cumsum(route_details$dist)
route_details$coord <- node_coords
route_details$id <- node_ids
print(route_details)

#here we select randomly parts from the route

gps_points<-route_details[sample(1:nrow(route_details), 10,replace=FALSE),] 

here i would like to change the weight of the graph parts based on selected gps coordinates. i get the way till the point of getting the gps coordinates but i just hang up mentally here to get back to the graph to change the weights there.

# Currently i can only Modify current route weight
E(gr)[r]$weight <- E(gr)[r]$weight * 2

Thank you for you help! Best regards.

Andreas
  • 397
  • 4
  • 18
  • 37

1 Answers1

1

The following script finds the ids of edges adjacent to a list of coordinates (wished.coord) so that you can modify the weights :

library(osmar)
library(igraph)
library(tidyr)
library(dplyr)

### Get data ----
src <- osmsource_api(url = "https://api.openstreetmap.org/api/0.6/")
muc_bbox <- center_bbox(11.575278, 48.137222, 1000, 1000)
muc <- get_osm(muc_bbox, src)

### Reduce to highways: ----
hways <- subset(muc, way_ids = find(muc, way(tags(k == "highway"))))
hways <- find(hways, way(tags(k == "name")))
hways <- find_down(muc, way(hways))
hways <- subset(muc, ids = hways)

#### Plot data ----
## Plot complete data and highways on top:
plot(muc)
plot_ways(muc, col = "lightgrey")
plot_ways(hways, col = "coral", add = TRUE)

### Define route start and end nodes: ----
id<-find(muc, node(tags(v %agrep% "Sendlinger Tor")))[1]
hway_start_node <-find_nearest_node(muc, id, way(tags(k == "highway"))) 
hway_start <- subset(muc, node(hway_start_node))

id <- find(muc, node(attrs(lon > 11.58 & lat > 48.15)))[1]
hway_end_node <- find_nearest_node(muc, id, way(tags(k == "highway")))
hway_end <- subset(muc, node(hway_end_node))

## Add the route start and and nodes to the plot:
plot_nodes(hway_start, add = TRUE, col = "red", pch = 19, cex = 2)
plot_nodes(hway_end, add = TRUE, col = "red", pch = 19, cex = 2)

### Create street graph ----
gr <- as.undirected(as_igraph(hways))

### Compute shortest route: ----
# Calculate route
route <- function(start_node,end_node) {
  get.shortest.paths(gr,
                     from = as.character(start_node),
                     to = as.character(end_node), 
                     mode = "all")[[1]][[1]]}
# Plot route
plot.route <- function(r,color) {
  r.nodes.names <- as.numeric(V(gr)[r]$name)
  r.ways <- subset(hways, ids = osmar::find_up(hways, node(r.nodes.names)))
  plot_ways(r.ways, add = TRUE, col = color, lwd = 2)
}
nways <-  1
numway <- 1
r <- route(hway_start_node,hway_end_node)

# Plot route

color <- colorRampPalette(c("springgreen","royalblue"))(nways)[numway]
plot.route(r,color)


## Route details ----
# Construct a new osmar object containing only elements 
# related to the nodes defining the route:
route_nodes <- as.numeric(V(gr)[r]$name)
route_ids <- find_up(hways, node(route_nodes))

osmar.route <- subset(hways, ids = route_ids)
osmar.nodes.ids <- osmar.route$nodes$attrs$id

# Extract the nodes’ coordinates,
osmar.nodes.coords <- osmar.route$nodes$attrs[, c("lon", "lat")]
osmar.nodes <- cbind(data.frame(ids = osmar.nodes.ids),
                     data.frame(ids_igraph = as.numeric(V(gr)[r]) ),
                     osmar.nodes.coords) 


## Find edges ids containing points of interest ----
wished.coords <- data.frame(wlon = c(11.57631),
                            wlat = c(48.14016)) 


# Calculate all distances
distances <- crossing(osmar.nodes,wished.coords) %>%
             mutate(dist = geosphere::distHaversine(cbind(lon,lat),cbind(wlon,wlat)))


# Select nodes below maximum distance :
mindist <- 50 #m

wished.nodes <- distances %>% filter(dist < mindist)

# Select edges incident to these nodes :
selected.edges <- unlist(incident_edges(gr,V(gr)[wished.nodes$ids_igraph]))

# Weight of selected edges
E(gr)[selected.edges]$weight
Waldi
  • 39,242
  • 6
  • 30
  • 78
  • thank you: i get an error Fehler in membership(communities) : Cannot calculate community membership – Andreas Jun 16 '20 at 21:44
  • on which command of the script do you get this error? – Waldi Jun 16 '20 at 21:46
  • # Calculate all distances distances <- crossing(osmar.nodes,wished.coords) %>% mutate(dist = geosphere::distHaversine(cbind(lon,lat),cbind(wlon,wlat))) – Andreas Jun 16 '20 at 21:47
  • did you load packages tidyr, is geosphere installed on your computer? – Waldi Jun 16 '20 at 21:52
  • Those are packages i load: library(osmar) library(igraph) library(tidyr) library(dplyr) library(geosphere) – Andreas Jun 16 '20 at 21:54
  • where are values for lon, lat in the line : distHaversine(cbind(lon,lat) are commming from? maybe this is the error? – Andreas Jun 16 '20 at 21:56
  • from osmar.nodes.coords. I propose to test a complete restart of RStudio to clean-up session – Waldi Jun 16 '20 at 21:59
  • I sourced script after restart : no error. Result = E(gr)[selected.edges]$weight [1] 9.138732 11.458915 14.286534 11.458915 12.379732 20.426907 17.470794 6.378125 17.470794 6.669311 3.180464 [12] 18.319578 9.897568 3.180464 9.138732 9.897568 6.669311 8.678559 18.443744 14.286534 6.378125 – Waldi Jun 16 '20 at 22:02
  • did a restart, the script works now. Thank you very much again Waldi. – Andreas Jun 16 '20 at 22:03
  • My pleasure @Andreas! – Waldi Jun 16 '20 at 22:05
  • this operation selected.edges <- unlist(incident_edges(gr,V(gr)[wished.nodes$ids_igraph])) if i have a lot of gps points is quite slow. Is there a way to optimize it? – Andreas Jun 17 '20 at 08:46
  • to be precise: this operation incident_edges(gr,V(gr)[wished.nodes$ids_igraph]) takes quite a long time. is there a way to optimize it? my wishlist has around 300 rows – Andreas Jun 17 '20 at 09:09
  • I tried to microbenchmark it on one node : microbenchmark::microbenchmark(incident_edges(gr,V(gr)[wished.nodes$ids_igraph[1]])) and get ~10ms per node, which would suggest 3s for 300 nodes. How long does it take? – Waldi Jun 17 '20 at 10:59
  • expr min lq mean median selected.edges <- unlist(incident_edges(gr_muc, V(gr_muc)[wished.nodes$ids_igraph[1]])) 50.62764 52.75964 60.83878 53.77717 uq max neval 58.8037 124.7251 100 – Andreas Jun 17 '20 at 11:20
  • nrow(wished.nodes) [1] 780 – Andreas Jun 17 '20 at 11:21
  • i stopped the time again: it took 46.763 sec or nrow(wished.nodes)=780 and gsize(gr)=85401 for incident_edges operation. Would be great if it could be optimized. – Andreas Jun 17 '20 at 12:42
  • Let us [continue this discussion in chat](https://chat.stackoverflow.com/rooms/216145/discussion-between-andreas-and-waldi). – Andreas Jun 17 '20 at 15:12
  • i came back to this solution because i think the following code :mindist <- 50 #m wished.nodes <- distances %>% filter(dist < mindist) does contain a logical error. So we have a dataframe of nodes and calculate the distances between the n-node and n+1 node. Then if the distance is smaller than min distance, we filter the node n+1 out and move to next node for compare after the deleted one and pick next compare. But the problem as i see it we should still check the distance between the first and the 3d node and so on till we find one that satisfies the condition and then move on? – Andreas Jul 22 '20 at 11:44
  • @Andreas, the code finds all the nodes on the path that are near (less than 50m) points of interest, and modifies the weight of the edges bound to these nodes. Crossing does a cartesian product : we compare all nodes on the path with all nodes available. I don't think there is a n / n+1 logic – Waldi Jul 22 '20 at 14:26
  • Waldi you a correct i think i mixed up my "problems" sorry for that. Currently im trying to thin out the GPS trace and i used this approach which lead me to beleve it is possible to use this same approach to filter the GPS trace.. which lead to wrong connclusion that there is a n+1 logic. Im realy sorry. This is where im currently stuck .. it looked similar : https://stackoverflow.com/questions/63034689/write-an-efficient-loop-to-compare-gps-coordinates – Andreas Jul 22 '20 at 14:29
  • OK! did you manage to create the 1km grid with the comment of @Limey, or are you still interested in a solution? – Waldi Jul 22 '20 at 14:39
  • Im absolutely interrested in a solution! Currently i just use a "hack" longitudes <- c(8.9771580802, 13.8350427083) latitudes <- c(47.2703623267, 50.5644529365) bounding_box <- matrix(c(longitudes, latitudes), nrow = 2, byrow = TRUE, dimnames = list(NULL, c("min", "max"))) projection <- "+proj=longlat" sp_box<-Spatial(bbox = bounding_box, proj4string = CRS(projection)) points_bavaria <- spsample(sp_box, 705000, type="regular") – Andreas Jul 22 '20 at 14:50