Estrazione di ristoranti in un raggio attorno ad un punto

script.R

library(osmdata)
library(sf)

# Punto a caso in Israele
coords <- c(34.8, 32)
coords <- st_sfc(st_point(x = coords), crs = 4326)
# Trasformazione nella proiezione Israeli TM Grid
# https://spatialreference.org/ref/epsg/israel-israeli-tm-grid/
coords_proj <- st_transform(coords, 2039)
# Ho proiettato così da creare un vero cerchio con un raggio in metri
circle_proj <- st_buffer(coords_proj, dist=10000) # 10000 m = 10 km di raggio
# Osmdata vuole che gli si parli in lat-lon quindi riproietto il cerchio in lat-lon
circle <- st_transform(circle_proj, 4326)

# Estrazione della bounding box
bb <- st_bbox(circle, crs = 4326)

# Query sul server ufficiale di OSM
q <- opq(as.vector(bb)) %>% # caveat! credo non gli piaccia la classe bbox e per questo faccio il cast in vector
  add_osm_feature("amenity", "restaurant") %>% 
  osmdata_sf()

# Ottengo punti e poligoni. Mi tengo solo i punti per comodità.
# Ora faccio l'intersezione con il mio cerchio iniziale
mypts <- st_intersection(q$osm_points, circle)

# Li disegno sulla mappa
library(leaflet)

leaflet() %>% addTiles() %>%
  addPolygons(data=circle) %>% 
  addMarkers(data=mypts) %>% 
  addScaleBar(options = scaleBarOptions(imperial=FALSE))