R 可视化：地图上的合理排斥点（蜂群？）

Question

我正在尝试复制大约这样的地图。

它描绘了分布在一个区域内的少量项目（学校）。对于输入，我有一个区域地图，每个区域都有一个数字。我想把它放在该地区周围的许多点上。如果它们不跨越区域边界扩散会更好，但简单地分布它们就足够了。区域内一些不错的排斥点可能会起作用。

Beeswarm plots 做了一些非常相似的事情，这可以在地图上完成吗？额外的问题 - 事实上我一直在寻找动画，但只能想出非常复杂的方法来做到这一点，以便随着时间的增加和 nrs 添加新点。

下面的代码将点放在地图上的质心中，并将数字作为大小。 （我无法将地图正确导出为单个文件，所以坐标有点乱，但原理是一样的。）

places = st_read("https://gist.githubusercontent.com/peeter-t2/9646a4169e993948fa97f6f503a0688b/raw/cb4e910bf153e51e3727dc9d1c73dd9ef86d2556/kih1897m.geojson", stringsAsFactors = FALSE)

schools <- read_tsv("https://gist.github.com/peeter-t2/34467636b3c1017e89f33284d7907b42/raw/6ea7dd6c005ef8577b36f5e84338afcb6c76b707/school_nums.tsv")
schools_geo <- merge(places,schools,by.x="KIHELKOND",by.y="Kihelkond") #94 matches

p<- schools_geo %>% 
  ggplot()+
  geom_sf(data=schools_geo)+
  geom_sf(data=st_centroid(schools_geo),aes(size=value))+
  theme_bw()
p

谢谢！

Answer 1

正如我在 cmets 中所指出的，当我在文件中读取时，它将 crs 设置为纬度/经度 (epsg: 4326)，而几何列是不同的 crs。我猜到正确的crs 是espg: 3301 并在此基础上继续进行，这似乎工作正常。

st_crs(schools_geo) <- 3301

我们可以使用 st_sample 来获取多边形内与我们的“值”列相关的点样本：

# we can set type = 'hexagonal', 'regular' or 'random'
school_pts <- schools_geo %>% st_sample(size = .$value, type = 'hexagonal')


schools_geo %>% 
  ggplot()+
  geom_sf()+
  geom_sf(data=school_pts, size = .8)+
  theme_bw()

这会产生下面的图，我认为它看起来很乱，因为st_sample 将点扩展到多边形的范围内。

像您发布的示例中那样，使每个多边形中的点更居中可能看起来会更好。为此，我们可以根据要在其中绘制的点数重新缩放多边形。在下面的代码中，如果多边形内部的点数最少 (1)，我将多边形缩小 90%，如果它们的点数最多 (27)，则缩小 20%。

# put values on scale between 0 and 1
scale_fact <- (max(schools_geo$value) -  schools_geo$value) / (max(schools_geo$value) - min(schools_geo$value)) 
# re-scale between 0.2 and 0.9
scale_fact <- scale_fact * (0.9 - 0.2) + 0.2
# reverse the scale 
scale_fact <-  max(scale_fact) + min(scale_fact) - scale_fact 

# apply the scale factor
schools_centroid <- st_geometry(st_centroid(schools_geo))
schools_geo_rescaled <- (st_geometry(schools_geo) - schools_centroid) * scale_fact + schools_centroid

school_pts <- schools_geo_rescaled %>% 
  st_sf(crs = 3301) %>% 
  bind_cols(value = schools_geo$value) %>%
  st_sample(size = .$value, type = 'hexagonal')


# plot
schools_geo %>% 
  ggplot()+
  geom_sf()+
  geom_sf(data=school_pts, size = .8)+
  theme_bw()

Answer 2

这不是一个简单的问题。我决定简化它，只选择一个区域而不是全部。理论上，该解决方案适用于您的所有领域。

我们首先导入我们的库

library(rgdal)
library(sf)
library(readr)
library(ggplot2)

我们使用建议的数据：

places <- st_read("https://gist.githubusercontent.com/peeter-t2/9646a4169e993948fa97f6f503a0688b/raw/cb4e910bf153e51e3727dc9d1c73dd9ef86d2556/kih1897m.geojson", stringsAsFactors = FALSE)

schools <- read_tsv("https://gist.github.com/peeter-t2/34467636b3c1017e89f33284d7907b42/raw/6ea7dd6c005ef8577b36f5e84338afcb6c76b707/school_nums.tsv")
schools_geo <- merge(places,schools,by.x="KIHELKOND",by.y="Kihelkond") #94 matches

我们选择一个状态

one <- places$geometry[[1]]

借助网格，我们将多边形分割成多个子多边形

grid <- st_make_grid(one, n = c(10, 10))
area <- st_area(grid)
grid <- st_as_sf(data.frame(ID=c(1:length(area)),
                            area=area,
                            geometry=grid))
tmp <- st_intersection(grid, one)
tmp$area <- st_area(tmp)

我们显示由小正方形组成的网格的所有质心

plot(st_geometry(tmp['area']))
plot(st_geometry(st_centroid(tmp['area'])),
     pch = 16, col = 'red', add = TRUE)

最后，我们只想保留您想要的分数，这在您的示例中相当于价值（学校数量）。

nbr <- 20
plot(st_geometry(one))
plot(st_geometry(st_centroid(tmp[order(tmp$area, decreasing = T),][1:nbr,])),
     pch = 16, col = 'red', add = TRUE)

希望对你有所帮助。