gis - Reverse cluster analysis; identifying empty space or a lack of density in R with longitude and latitude? -

i'm working on project have large amount of points , looking identify regions (defined lack of clustering) density of these points statistically less relative others. visual enough have many points difficult tell these empty spaces , density heat map doesn't me 0 in on smaller regions. maybe i'm missing simpler here, hoping can @ least send me in right direction of look. below reproducible sample quick , dirty lets take these points open data , map them borough file nyc:

#libraries--------------------------  library(ggplot2) library(ggmap) library(sp) library(jsonlite) library(rjsonio) library(rgdal)  #call api data--------------------------  df =  fromjson("https://data.cityofnewyork.us/resource/24t3-xqyv.json?$query= select lat, long_") df <- data.frame(t(matrix(unlist(df),nrow=length(unlist(df[1]))))) names(df)[names(df) == 'x2'] = 'x' names(df)[names(df) == 'x1'] = 'y' df = df[, c("x", "y")] df$x = as.numeric(as.character(df$x)) df$y = as.numeric(as.character(df$y)) df$x = round(df$x, 4) df$y = round(df$y, 4) df$x[df$x < -74.2] = na df$y[df$y < 40.5] = na df = na.omit(df)  #map data----------------------------   cd = readogr("nybb.shp", layer = "nybb") cd = sptransform(cd, crs("+proj=longlat +datum=wgs84")) cd_f = fortify(cd)   #map data nyc = ggplot() +   geom_polygon(aes(x=long,                     y=lat, group=group), fill='grey',                 size=.2,color='black', data=cd_f, alpha=1)   nyc + geom_point(aes(x = x, y = y), data = df, size = 1)  #how go finding empty spaces? regions there no clusters? 

in case there aren't lot of points sake of demonstration, how i:

1. identify pockets of low density
2. potentially draw polygon boundaries on pockets?

appreciate help!

one way polygonal areas of low density construct dirichlet/voronoi tesselation , choose largest ones.

below use spatstat , deldir (loaded spatstat) this. not fast many more points don't know how go.

to use results in ggmap , other spatial packages can convert owin , ppp spatial classes sp , use sptransform lat, long coordinates.

first load packages:

library(maptools) library(spatstat) library(jsonlite) 

map , points in coordinates of shapefile (note read in data local files downloaded www):

cd = readogr("nybb.shp", layer = "nybb") #> ogr data source driver: esri shapefile  #> source: "nybb.shp", layer: "nybb" #> 5 features #> has 4 fields df <- fromjson("nyc_data.json") df <- as.data.frame(matrix(as.numeric(unlist(df)), ncol = 2, byrow = true)) df <- df[, c(2, 1)] names(df) <- c("x", "y") df <- df[df$x > -74.2 & df$y > 40.5, ] coordinates(df) <- ~x+y proj4string(df) <- crs("+proj=longlat +datum=wgs84") df2 <- sptransform(df, proj4string(cd)) 

switch spatstat classes:

cd2 <- as(cd, "spatialpolygons") w <- as(cd2, "owin")  x <- as(df2, "ppp") window(x) <- w  plot(x, main = "") 

compute dirichlet tessellation , areas , plot tessellation:

d <- dirichlet(x) #> warning: 96 duplicated points removed <- tile.areas(d) plot(d, main = "") 

combine n_areas biggest areas of tessellation:

n_areas <- 30 empty <- tess(tiles = d$tiles[tail(order(a), n = n_areas)]) empty2 <- as.owin(empty) 

plot result:

plot(w, main = "") plot(empty2, col = "red", add = true) plot(x, add = true)