Анализ обратного кластера; идентифицируя пустое пространство или отсутствие плотности в R с долготой и широтой?

Question

Я работаю над проектом, где у меня очень много точек, и я ищу для определения областей (определяемых отсутствием кластеризации), где плотность этих точек статистически значительно меньше по сравнению с другими. Обычно визуального было бы достаточно, но у меня так много точек, что трудно сказать, где эти пустые пространства и плотность тепловой карты плотности не помогает мне нуля в меньших областях. Может быть, мне здесь что-то не так просто, но я надеюсь, что кто-то может по крайней мере отправить меня в правильном направлении, где искать. Ниже воспроизводимый образец быстро и грязно позволяет принимать эти точки из открытых данных и отобразить их в файл р-н для Нью-Йорка:

#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 would I go about finding the empty spaces? That is the regions where there are no clusters? 

В этом случае есть не много очков, но ради демонстрации, как бы я:

1. определить очаги низкой плотности
2. потенциально начертить границы полигонов на этих карманах?

Цените помощь!

Answer 1

Одним из способов получения многоугольных участков с низкой плотностью было бы построение тесселяции Dirichlet/Voronoi и выбор наиболее крупных.

Ниже я использую spatstat и deldir (загружен spatstat) для этого. Это не так быстро, поэтому со многими другими точками я не знаю, как хорошо это будет go.

Чтобы использовать результаты в ggmap и других пространственных пакетов вы можете конвертировать назад от owin и ppp к пространственным классам от sp и использовать spTransform получить лат, длинные координаты.

сначала загрузить пакеты:

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

Карты и точки в координатах шейпа (обратите внимание, я прочитал в данных из локальных файлов, загруженных из WWW):

cd = readOGR("nybb.shp", layer = "nybb") 
#> OGR data source with driver: ESRI Shapefile 
#> Source: "nybb.shp", layer: "nybb" 
#> with 5 features 
#> It 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)) 

Переключить на spatstat классов:

cd2 <- as(cd, "SpatialPolygons") 
W <- as(cd2, "owin") 

X <- as(df2, "ppp") 
Window(X) <- W 

plot(X, main = "") 

тесселяции и области

Compute Дирихле и сюжет тесселяции:

d <- dirichlet(X) 
#> Warning: 96 duplicated points were removed 
a <- tile.areas(d) 
plot(d, main = "") 

Объедините n_areas большие участки тесселяции:

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

земля Результат:

plot(W, main = "") 
plot(empty2, col = "red", add = TRUE) 
plot(X, add = TRUE)