R中的二进制密度图

Binary density map in R

我有一个大型数据集,其中包含家庭位置的 x 和 y 坐标对。

我想要的是创建一个区分密集区域和不密集区域的二进制密度图(确切的阈值将通过启发式定义)。因此,我想说如果一个坐标对在某个区域内,它就被认为是密集的。在这种情况下,我不确定克里金法或反距离加权是否是合适的方法。

到目前为止,我的方法如下:

library(ggplot2)

m <- ggplot(df, aes(x = df$V1, y = df$V2)) +
  geom_point() 
m + geom_density_2d()

不幸的是,这并不能真正让我提取 binary map,而只能显示密度。一个可以想到的解决方案是第三列,它将坐标对分为密集和非密集类别。

非常感谢您的帮助!

数据 数据摘录(总计 162'000 行):

df <- structure(list(V1 = c(2690503, 2689797, 2690685, 2690685, 2689409, 
2689409, 2689141, 2689141, 2689141, 2689141, 2690515, 2690515, 
2689474, 2689148, 2689148, 2689148, 2689148, 2689148, 2689148, 
2690199, 2689760, 2689473, 2689198, 2689198, 2689558, 2690020, 
2690020, 2689612, 2689132, 2689208, 2689208, 2690244, 2689614, 
2689614, 2690264, 2690264, 2689625, 2690315, 2689720, 2689720, 
2689720, 2691037, 2691037, 2691037, 2691037, 2689433, 2689433, 
2690715, 2690715, 2690715, 2689560, 2689655, 2689563, 2690240, 
2690240, 2678781, 2689498, 2689206, 2690521, 2689662, 2689662, 
2690027, 2690027, 2689383, 2689383, 2690588, 2688999, 2689397, 
2689725, 2689725, 2689100, 2689100, 2689100, 2689100, 2689906, 
2690529, 2690529, 2690199, 2690164, 2689638, 2689638, 2689498, 
2689093, 2689093, 2690502, 2689740, 2689153, 2690536, 2689027, 
2689950, 2690959, 2690959, 2690486, 2690486, 2689794, 2689307, 
2690010, 2690010, 2689599, 2689599, 2689599, 2689464, 2689464, 
2689464, 2689464, 2689711, 2689711, 2690222, 2690602, 2690602, 
2689790, 2689790, 2689404, 2689428, 2689428, 2689387, 2688960, 
2688960, 2689258, 2689258, 2689258, 2689355, 2690795, 2689521, 
2689954, 2689954, 2688926, 2689543, 2689363, 2689363, 2689186, 
2689186, 2689571, 2689571, 2689970, 2689970, 2689675, 2689498, 
2690941, 2690941, 2689060, 2689442, 2690122, 2690725, 2690725, 
2689419, 2689366, 2689366, 2689097, 2689332, 2690123, 2690123, 
2690994, 2690994, 2690180, 2690180, 2689706, 2689706, 2689612, 
2690074, 2688961, 2688961, 2689692, 2689260, 2689419, 2689419, 
2689146, 2690083, 2690625, 2690167, 2690167, 2689540, 2689540, 
2689512, 2689512, 2690469, 2689720, 2689711, 2690874, 2690072, 
2690072, 2690072, 2688946, 2689502, 2689431, 2689531, 2689131, 
2689131, 2690257, 2690001, 2689608, 2689843, 2689502, 2689773, 
2689773, 2689507, 2690060, 2678781, 2689500, 2689260), V2 = c(1254816, 
1254916, 1254061, 1254061, 1255542, 1255542, 1255220, 1255220, 
1255220, 1255220, 1254872, 1254872, 1255561, 1255199, 1255199, 
1255199, 1255199, 1255199, 1255199, 1255390, 1255667, 1255233, 
1255830, 1255830, 1255029, 1254812, 1254812, 1255297, 1255391, 
1255728, 1255728, 1254961, 1255385, 1255385, 1255149, 1255149, 
1255704, 1255312, 1254949, 1254949, 1254949, 1253836, 1253836, 
1253836, 1253836, 1255130, 1255130, 1253886, 1253886, 1253886, 
1255124, 1254928, 1255858, 1255267, 1255267, 1237314, 1255231, 
1255426, 1254796, 1255315, 1255315, 1255231, 1255231, 1255065, 
1255065, 1254882, 1255504, 1255493, 1255279, 1255279, 1256005, 
1256005, 1256005, 1256005, 1255418, 1254909, 1254909, 1255390, 
1255233, 1255716, 1255716, 1255231, 1255787, 1255787, 1253745, 
1255672, 1255827, 1254775, 1255813, 1255187, 1254105, 1254105, 
1255155, 1255155, 1255128, 1255623, 1255448, 1255448, 1255397, 
1255397, 1255397, 1255353, 1255353, 1255353, 1255353, 1255306, 
1255306, 1254824, 1254771, 1254771, 1255170, 1255170, 1255380, 
1255919, 1255919, 1255204, 1255885, 1255885, 1256001, 1256001, 
1256001, 1255552, 1254091, 1255334, 1255052, 1255052, 1255609, 
1254960, 1255090, 1255090, 1255426, 1255426, 1255140, 1255140, 
1254886, 1254886, 1255464, 1255231, 1254052, 1254052, 1255454, 
1255219, 1255265, 1254818, 1254818, 1255145, 1255407, 1255407, 
1255180, 1255187, 1255179, 1255179, 1253875, 1253875, 1255001, 
1255001, 1255573, 1255573, 1255232, 1255023, 1255458, 1255458, 
1255305, 1255453, 1255145, 1255145, 1255275, 1255209, 1254931, 
1255447, 1255447, 1255003, 1255003, 1255760, 1255760, 1254749, 
1254949, 1255306, 1253851, 1255286, 1255286, 1255286, 1255599, 
1255183, 1255580, 1255320, 1255460, 1255460, 1254804, 1255360, 
1255731, 1255470, 1255183, 1255131, 1255131, 1255081, 1255245, 
1237314, 1255586, 1255419)), row.names = c(NA, 200L), class = "data.frame")

一个类似的问题has been asked in this thread。我正在修改给出的答案。

我正在使用 MASS:kde2d 制作您实际上也提到的数据框。请注意,此函数将给出矩阵 z ,它 - 像往常一样 - 按列排列。因此,您需要将 y 坐标重复 each,将 x 重复 times

library(tidyverse)

# Manually estimate the density for each coordinate value with MASS::kde2d
# I scale the data first and remove this weird outlier in your data 
scaled_df <- df %>% scale() %>% as.data.frame() %>% filter (V1 > -7)

kde2d_est <- MASS::kde2d(scaled_df$V1, scaled_df$V2)

res <- data.frame(x = rep(kde2d_est$x, times = dim(kde2d_est$z)[1]),
                  y = rep(kde2d_est$y, each = dim(kde2d_est$z)[2]),
                  density = as.numeric(kde2d_est$z))

ggplot() +
  geom_point(data = scaled_df, aes(V1, V2)) +
  geom_tile(data = res, aes(x = x, y = y, fill = density), alpha = 0.5)

看起来已经不错了

现在让我们制作 "binary density",任意截断。您不需要像我链接到的答案中那样 'cuts' 。只需使用逻辑语句即可。

arbitrary_cut <- 1
res_bin <- res %>% mutate(binary_dens = density < arbitrary_cut)

ggplot() +
  geom_point(data = scaled_df, aes(V1, V2)) +
  geom_tile(data = res_bin, aes(x = x, y = y, fill = binary_dens), alpha = 0.5)

reprex package (v0.3.0)

于 2020-04-05 创建