mavat_to_net.R

library(sf)
library(tidyverse)
library(sfnetworks)
library(tidygraph)
library(lwgeom)
library(mapview)
library(igraph)
library(qgisprocess)
library(nngeo)
library(dbscan)
library(leafsync)
# options(qgisprocess.path = "C:/OSGeo4W64/bin/python-qgis.bat")
# qgis_configure()
st_layers("C:/Users/idshk/OneDrive/idos_shit/kollelani.gpkg")
rehovot <- st_read("C:/Users/idshk/OneDrive/idos_shit/kollelani.gpkg",layer = "rehovot_east",crs=2039)
roads <- rehovot %>% filter(str_detect(MAVAT_NAME,"דרך"))
roads_valid <- st_make_valid(roads)

roads_buf <- st_buffer(roads_valid,1)
roads_dissolve <- summarise(roads_buf)
roads_debuf <- st_buffer(roads_dissolve,-1)
roads_single <- st_cast(roads_debuf,"POLYGON")
rodas_lines <- roads_single[1,] %>% st_cast("LINESTRING") 
roads_densify <- st_line_sample(rodas_lines,density = 1/10)  
roads_points <- st_cast(roads_densify,"POINT") 
bbox <- st_as_sfc(st_bbox(st_buffer(st_cast(roads_points,"MULTIPOINT"),100)))
voronoi <- st_collection_extract(st_voronoi(st_union(roads_points), bbox))
lined_voronoi <- st_cast(voronoi,"LINESTRING") 
segment_voronoi <- st_segments(lined_voronoi)
to_keep <- t(st_contains(roads_dissolve,segment_voronoi,sparse = F))
streets <- st_sf(segment_voronoi)[to_keep,] 
to_remove <- streets %>% 
  mutate(start = st_startpoint(segment_voronoi),
         end = st_endpoint(segment_voronoi),
         vec = map2(start,end,~sort(c(st_coordinates(.x),st_coordinates(.y))))) %>% 
  pull(vec) %>% 
  duplicated()

net <- as_sfnetwork(streets[!to_remove,],directed = F)
m1 <- net %>%activate(edges) %>%  st_as_sf() 
m2 <- net %>% st_as_sf()
m3 <- mapview(m1) + m2
smoothed = convert(net, to_spatial_smooth)
m4 <- smoothed %>%activate(edges) %>%  st_as_sf() %>% select(-3)
m5 <- smoothed %>% st_as_sf() 
m6 <- mapview(m4) + m5
node_coords = smoothed %>%
  activate("nodes") %>%
  st_coordinates()
clusters = dbscan(node_coords, eps = 35, minPts = 1)$cluster
clustered = smoothed %>%
  activate("nodes") %>%
  mutate(cls = clusters)
clustered = clustered %>%
  mutate(cmp = group_components())
contracted = convert(
  clustered,
  to_spatial_contracted,
  cls, cmp,
  simplify = TRUE
)
m7 <- contracted %>%activate(edges) %>%  st_as_sf() %>% select(-4)
m8 <- contracted %>% st_as_sf() %>% select(-1)
m9 <- mapview(m7) + m8
subdivision = convert(contracted, to_spatial_subdivision)
m10 <- subdivision %>%activate(edges) %>%  st_as_sf() %>% select(-4)
m11 <- subdivision %>% st_as_sf() %>% select(-1)
m12 <- mapview(m10) + m11
node_coords2 = subdivision %>%
  activate("nodes") %>%
  st_coordinates()
clusters2 = dbscan(node_coords2, eps = 35, minPts = 1)$cluster
clustered2 = subdivision %>%
  activate("nodes") %>%
  mutate(cls = clusters2)
clustered2 = clustered2 %>%
  mutate(cmp = group_components())
contracted2 = convert(
  clustered2,
  to_spatial_contracted,
  cls, cmp,
  simplify = TRUE
)
m13 <- contracted2 %>%activate(edges) %>%  st_as_sf() %>% select(-4)
m14 <- contracted2 %>% st_as_sf() %>% select(-1)
m15 <- mapview(m13) + m14
smoothed2 = convert(contracted2, to_spatial_smooth)
m16 <- smoothed2 %>%activate(edges) %>%  st_as_sf() %>% select(-3) %>% st_simplify(preserveTopology = TRUE,dTolerance = 35)
m17 <- smoothed2 %>% st_as_sf() %>% select(-1)
m18 <- mapview(m16) + m17
m18