new_buf_approach.R

library(tidyverse)
library(magrittr)
library(sf)
library(lwgeom)
library(sfnetworks)
library(osmdata)
library(ggspatial)
library(mapview)
library(tidygraph)
library(SpatialPosition)
library(parallel)
library(concaveman)
library(ggspatial)
library(TSP)
library(igraph)
# download data
q1 <- opq(c(35.18,31.75,35.22,31.77)) %>%
  add_osm_feature(key = 'highway') %>% 
  osmdata_sf()

# preprocess data
# roundabouts are coded as polygons
center_jlm <- bind_rows(st_transform(q1$osm_lines,2039),st_transform(q1$osm_polygons,2039)) %>% 
  st_cast("LINESTRING")
# split network in middlepoints
nodes_to_split_by <- center_jlm %>%
  st_cast("POINT") %>% 
  group_by(osm_id) %>% 
  filter(row_number() == 1| row_number() == max(row_number())) %>% 
  st_geometry() %>% 
  st_sf() %>% 
  distinct() %>% 
  st_snap(center_jlm,tolerance = 1e-9)
center_jlm_break <- st_split(center_jlm$geometry, nodes_to_split_by$geometry) %>% 
  st_collection_extract("LINESTRING") %>% 
  st_sf()
# create sfnetworks object
net <- as_sfnetwork(center_jlm_break, directed = FALSE)
net <- net %>% 
  activate(edges) %>% 
  mutate(len = as.numeric(st_length(geometry)))

net <- net %>% 
  activate(nodes) %>% 
  mutate(rn = row_number())

nodes <- net %>% 
  activate(nodes) %>% 
  st_as_sf()

# distances calculation
sparse_shortest_path_dist <- distances(net,weights = E(net)$len)
row.names(sparse_shortest_path_dist) <- 1:nrow(sparse_shortest_path_dist)
colnames(sparse_shortest_path_dist) <- row.names(sparse_shortest_path_dist)
sparse_shortest_path_dist <- sparse_shortest_path_dist %>% 
  as.data.frame() %>% 
  rownames_to_column("pointA") %>% 
  gather(pointB,shortest_path_dist,-pointA) 

sparse_dist <- CreateDistMatrix(nodes,nodes) %>% 
  as.data.frame() %>% 
  rownames_to_column("pointA") %>% 
  gather(pointB,dist,-pointA) 

# pseudo-isodist - should use rgeos::gInterpolate() to get more accurate distances on lines
buf <- 1000
example_node <-1229
dist_final <- sparse_dist %>% 
  left_join(sparse_shortest_path_dist, by = c("pointA","pointB")) %>% 
  filter(shortest_path_dist < buf) %>% 
  mutate(pointA = pointA %>%  as.numeric(),
         pointB = pointB %>% as.numeric())


dist_final %>% 
  filter(pointA == example_node) %>% 
  left_join(nodes,by= c("pointA"="rn")) %>% 
  left_join(nodes,by= c("pointB"="rn")) %>% 
  st_sf() %>% 
  group_by(pointA) %>% 
  summarise(convex = st_combine(geometry.y)) %>% 
  mutate(concave = st_zm(concaveman(.)$polygons)) %>% 
  st_convex_hull() %>% 
  as_tibble() %>% 
  left_join(nodes,by= c("pointA"="rn")) %>% 
  mutate(buffer = st_buffer(geometry,buf)) %>% 
  gather(type,geom,-pointA) %>%
  st_sf(crs=2039) %>% 
  ggplot(aes(color = type)) + 
  annotation_map_tile(zoom = 15)+
  geom_sf(fill = NA,size = 2) + 
  scale_size(guide="none")

# tsp
set.seed(12)
samp <- sample(dist_final$pointA,10)

stops <- sparse_dist %>% 
  left_join(sparse_shortest_path_dist, by = c("pointA","pointB")) %>% 
  mutate(pointA = pointA %>%  as.numeric(),
         pointB = pointB %>% as.numeric()) %>% 
  filter(pointA %in% samp, pointB %in% samp) %>% 
  select(-dist) %>% 
  spread(pointB,shortest_path_dist) %>% 
  column_to_rownames("pointA") %>% 
  as.matrix() %>% 
  as.TSP()

stops <- insert_dummy(stops,label = "cut")

sol <- solve_TSP(stops)

res <- matrix(names(cut_tour(sol, "cut")),ncol = 1) %>% 
  as_tibble() %>% 
  mutate(V1 = as.numeric(V1),
         to = lead(V1)) %>% 
  rename(from = V1) %>% 
  slice(-nrow(.)) %>% 
  mutate(path = map2(from, to, ~st_shortest_paths(net,.x,.y)))

plot_tsp <- lapply(res$path,function(x){
  x$vpath%>% 
    unlist() %>% 
    as.matrix(ncol=1) %>% 
    as_tibble() %>% 
    mutate(to = lead(V1)) %>% 
    slice(-nrow(.))
}) %>% 
  bind_rows() %>% 
  mutate(from = unlist(map2(V1,to,~min(.x,.y))),
         to = unlist(map2(V1,to,~max(.x,.y)))) %>% 
  select(from,to) %>% 
  left_join(net %>% activate(edges) %>% st_as_sf()) %>% 
  st_sf() 
tsp_points_order <- net %>% 
  activate(nodes) %>% 
  as_tibble() %>% 
  filter(rn %in% as.numeric(names(cut_tour(sol, "cut")))) %>% 
  arrange(factor(rn, levels = as.numeric(names(cut_tour(sol, "cut"))))) %>% 
  mutate(stop_num = row_number())
plot_tsp %>% 
  ggplot() + 
  annotation_map_tile(zoom = 15)+
  geom_sf() +
  geom_sf_label(data = tsp_points_order,mapping =aes(label = stop_num))