此外,这是一个基于 sf 和 tidyverse 的解决方案:
使用 sf,您可以使用以下命令制作一个正方形网格st_make_grid()
功能。在这里,我将在圣何塞的边界框上创建一个 2 公里长的网格,然后将其与圣何塞的边界相交。请注意,我正在投影到 UTM 区域 10N,因此我可以指定以米为单位的网格大小。
library(tigris)
library(tidyverse)
library(sf)
options(tigris_class = "sf", tigris_use_cache = TRUE)
set.seed(1234)
sj <- places("CA", cb = TRUE) %>%
filter(NAME == "San Jose") %>%
st_transform(26910)
g <- sj %>%
st_make_grid(cellsize = 2000) %>%
st_intersection(sj) %>%
st_cast("MULTIPOLYGON") %>%
st_sf() %>%
mutate(id = row_number())
接下来,我们可以生成一些随机犯罪数据st_sample()
并绘制它来看看我们正在处理什么。
thefts <- st_sample(sj, size = 500) %>%
st_sf()
assaults <- st_sample(sj, size = 200) %>%
st_sf()
plot(g$geometry)
plot(thefts, add = TRUE, col = "red")
然后可以将犯罪数据连接到空间网格中st_join()
。我们可以绘图来检查我们的结果。
theft_grid <- g %>%
st_join(thefts) %>%
group_by(id) %>%
summarize(num_thefts = n())
plot(theft_grid["num_thefts"])
然后我们可以对攻击数据执行相同的操作,然后将两个数据集连接在一起以获得所需的结果。如果您有大量犯罪数据集,可以修改这些数据集以使其在某些变化范围内工作purrr::map()
.
assault_grid <- g %>%
st_join(assaults) %>%
group_by(id) %>%
summarize(num_assaults = n())
st_geometry(assault_grid) <- NULL
crime_data <- left_join(theft_grid, assault_grid, by = "id")
crime_data
Simple feature collection with 190 features and 3 fields
geometry type: GEOMETRY
dimension: XY
bbox: xmin: 584412 ymin: 4109499 xmax: 625213.2 ymax: 4147443
epsg (SRID): 26910
proj4string: +proj=utm +zone=10 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
# A tibble: 190 x 4
id num_thefts num_assaults geometry
<int> <int> <int> <GEOMETRY [m]>
1 1 2 1 POLYGON ((607150.3 4111499, 608412 4111499, 608412 4109738,…
2 2 4 1 POLYGON ((608412 4109738, 608412 4111499, 609237.8 4111499,…
3 3 3 1 POLYGON ((608412 4113454, 608412 4111499, 607150.3 4111499,…
4 4 2 2 POLYGON ((609237.8 4111499, 608412 4111499, 608412 4113454,…
5 5 1 1 MULTIPOLYGON (((610412 4112522, 610412 4112804, 610597 4112…
6 6 1 1 POLYGON ((616205.4 4113499, 616412 4113499, 616412 4113309,…
7 7 1 1 MULTIPOLYGON (((617467.1 4113499, 618107.9 4113499, 617697.…
8 8 2 1 POLYGON ((605206.8 4115499, 606412 4115499, 606412 4114617,…
9 9 5 1 POLYGON ((606412 4114617, 606412 4115499, 608078.2 4115499,…
10 10 1 1 POLYGON ((609242.7 4115499, 610412 4115499, 610412 4113499,…
# ... with 180 more rows