反向聚类分析；用经度和纬度识别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 加载）来执行此操作。 它不是那么快，所以有更多的点我不知道它会有多好 去吧。

要使用ggmap 和其他空间包中的结果，您可以转换 从owin 和ppp 回到sp 的空间类并使用 spTransform 获取经纬度坐标。

首先加载包：

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

shapefile 坐标中的地图和点（注意我从 从 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 = "")

计算狄利克雷曲面细分和面积并绘制曲面细分：

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)