leaflet_thiessen.R

library(leaflet)
library(rgeos)
library(sp)
library(rgdal)
library(spatstat)
library(maptools)

leaflet_thiessen <- function(data, id, long, lat, epsg_code) {
  
  coords <- cbind(data[[long]], data[[lat]])
  
  ## Spatial points w/ the WGS84 datum
  
  projection <- paste0("+init=epsg:", as.character(epsg_code))
  
  sp_data <- SpatialPointsDataFrame(coords = coords, data = data, proj4string = CRS("+proj=longlat +datum=WGS84"))
  
  sp_data_proj <- spTransform(sp_data, CRS(projection))
  
  proj_coords <- sp_data_proj@coords
  
  window <- owin(range(proj_coords[,1]), range(proj_coords[,2]))
  
  d <- dirichlet(as.ppp(proj_coords, window))
  
  dsp <- as(d, "SpatialPolygons")
  
  dsp_df <- SpatialPolygonsDataFrame(dsp, 
                                     data = data.frame(id = 1:length(dsp@polygons)))
  
  proj4string(dsp_df) <- CRS(projection)
  
  dsp_df$area <- round((gArea(dsp_df, byid = TRUE) / 1000000), 1)
  
  dsp_xy <- spTransform(dsp_df, CRS("+proj=longlat +datum=WGS84"))
  
  
  leaflet() %>%
    addMarkers(data = data, 
               lat = data[[lat]], 
               lng = data[[long]], 
               popup = data[[id]]) %>%
    addPolygons(data = dsp_xy, 
                color = "green", 
                fill = "green", 
                popup = paste0("Area: ", 
                               as.character(dsp_xy$area), 
                               " square km")) %>%
    addTiles()
  
}