code to generate this map: https://rpubs.com/idshklein/jlm_walk_index_v2

buf_vs_cha.R

library(tidyverse)
library(magrittr)
library(sf)
library(lwgeom)
library(sfnetworks)
library(osmdata)
library(ggspatial)
library(mapview)
library(tidygraph)
q1 <- opq(c(35.18,31.75,35.22,31.77)) %>%
  add_osm_feature(key = 'highway') %>% 
  osmdata_sf()
# roundabouts are coded as polygons
center_jlm <- bind_rows(st_transform(q1$osm_lines,2039),st_transform(q1$osm_polygons,2039)) %>% 
  st_cast("LINESTRING")

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()

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())
number_of_rows <- net %>% 
  activate(nodes) %>% 
  as_tibble() %>% 
  nrow()
# improve algorithm - check shortest paths of only necessary neighbors
indices <- function(nrow,buf_dist){
  buf <- net %>% 
    activate(nodes) %>%
    st_as_sf() %>% 
    magrittr::extract(nrow,) %>% 
    st_buffer(buf_dist) %>% 
    mutate(area = as.numeric(st_area(geometry)),
           type = "buf")
  chp <- buf %>% 
    st_intersection(net) %>% 
    st_union() %>% 
    st_convex_hull() %>% 
    st_sf() %>% 
    mutate(rn = buf$rn,
           area = as.numeric(st_area(geometry)),
           type = "chp")
  
  cha <- net %>% 
    activate(nodes) %>% 
    st_as_sf() %>% 
    st_intersection(buf) %>% 
    mutate(dist = unlist(map2(rn,rn.1,~st_network_distance(net,.x,.y,weights = "len")[1]))) %>% 
    filter(dist < buf_dist) %>% 
    st_union() %>% 
    st_convex_hull() %>% 
    st_sf() %>% 
    mutate(rn = buf$rn,
           area = as.numeric(st_area(geometry)),
           type = "cha")
  print(paste0(nrow,"/",number_of_rows))
  bind_rows(buf,chp,cha) %>% 
    st_drop_geometry() %>% 
    spread(type,area)
}
indices_list <- vector("list", length = number_of_rows)
for(i in 1:number_of_rows){
  indices_list[[i]] <- indices(i,1000) 
}

df <- bind_rows(indices_list) %>% 
  left_join(net %>% 
              activate(nodes) %>% st_as_sf(), by = "rn") %>% 
  st_as_sf()
roundabouts <- center_jlm %>% filter(junction == "roundabout")
df %>% 
  mutate(cha_buf = cha/buf,
         chp_buf = chp/buf,
         cha_chp = cha/chp) %>%
  # filter(cha == 0) %>% 
  mapview(zcol = "cha_buf")  + center_jlm_break