There is an excellent blog post from a very curious and original person that made geographic vonoroni maps in R that I’ve adapted to use with OSM attributes and geometery.

Blog post: https://rstudio-pubs-static.s3.amazonaws.com/97256_125b33cf42414eda8797f53ac632628b.html

Code for doing this in OSM with a boundary shapefile.

library(sp)
library(rgdal)
library(deldir)
library(ggplot2)
library(ggthemes)
library(osmdata)
library(leaflet)
library(dplyr)
library(rgeos)
library(sf)
library(raster)

SPointsDF_to_voronoi_SPolysDF <- function(sp) {

  # tile.list extracts the polygon data from the deldir computation
  vor_desc <- tile.list(deldir(sp@coords[,1], sp@coords[,2]))

  lapply(1:(length(vor_desc)), function(i) {

    # tile.list gets us the points for the polygons but we
    # still have to close them, hence the need for the rbind
    tmp <- cbind(vor_desc[[i]]$x, vor_desc[[i]]$y)
    tmp <- rbind(tmp, tmp[1,])

    # now we can make the Polygon(s)
    Polygons(list(Polygon(tmp)), ID=i)

  }) -> vor_polygons

  # hopefully the caller passed in good metadata!
  sp_dat <- sp@data

  # this way the IDs _should_ match up w/the data & voronoi polys
  rownames(sp_dat) <- sapply(slot(SpatialPolygons(vor_polygons),
                                  'polygons'),
                             slot, 'ID')

  SpatialPolygonsDataFrame(SpatialPolygons(vor_polygons),
                           data=sp_dat)

}

### load Anchorage bourdary file
anc_boundary <- st_read(
  "Municipality_of_Anchorage_Boundary.shp")
anc_boundary <- as(anc_boundary, Class = "Spatial") # conver sf to sp

# query overpass API
x <- opq(bbox = c(-150.090972, 61.470444, -149.321814, 60.034781)) %>% 
  add_osm_feature(key = 'amenity', value = "place_of_worship", value_exact = FALSE) %>%
  osmdata_sf()

### Make vonroni polygons
pow_points <- st_centroid(x$osm_polygons)
pow_points <- as(pow_points, Class = "Spatial")
vor <- SPointsDF_to_voronoi_SPolysDF(pow_points)
vor <- intersect(anc_boundary, vor)
vor_df <- fortify(vor)

### Make Popup
y <- x$osm_polygons %>% 
  mutate(name = ifelse(!is.na(website), paste0('<a href="', website, '">', name, '</a>'), name),
         religion = paste0("religion: ", religion),
         denomination = ifelse(denomination != "none", paste0("denomination: ", denomination), ""),
         phone = ifelse(!is.na(website), paste0("phone #: ", phone), ""),
         email = ifelse(!is.na(website), paste0("email: ", email), "")) %>%
  mutate(for_popup = paste(name, phone, email, sep = "<br>"))

plot(intersect(anc_boundary, vor))


(map <- leaflet() %>% addTiles() %>% 
  setView(lng = -149.845, lat = 61.15, zoom = 12) %>%
  addPolygons(data = vor, fill = TRUE, fillOpacity = 0.1, weight = 1) %>% 
  addPolygons(data = y, color = "black", popup = ~for_popup, stroke = TRUE, fill = TRUE, weight = 3, 
              fillOpacity = 0.5, fillColor = "white",  options = popupOptions(autoPan = TRUE, closeButton = FALSE)))


library(htmlwidgets)

saveWidget(map, file = "anc_places_of_worship.html")