根据时间 window 计算多边形内的空间点
Counting spatial points within a polygon according to a temporal window
概览
使用 R，我想根据特定标准（时间window）计算多边形内的点数。
我有以下数据：
- 包含调查访谈日期的地理定位调查数据。因此，我能够准确指出每项调查的时间和地点，并将它们绘制在美国各地。
- 有关美国各地政治集会的地理定位数据。这些还包括日期。
使用 QGIS，我在每个调查受访者周围创建了一组 50 英里的圆形缓冲区。我的目标是计算在采访前的特定时间范围内每个受访者的“缓冲区”内的政治集会次数。 QGIS 中创建的 50 英里缓冲区保留了原始数据的所有变量，包括采访日期。
数据
我使用 QGIS 创建了一些 mock shapefiles 日期和位置以帮助复制。
方法
我正在尝试使用 GISTools::poly.counts 来计算不同时间 windows（30 天、90 天等）内的集会次数。
通常，要计算多边形内的点数，我会简单地使用：
count <- GISTools::poly.counts(rallies, buffer)
这为我提供了每个缓冲区内发生的集会总数，但不允许我指定时间 windows。例如，如果能够统计调查访谈前 30 天和访谈前 90 天缓冲区内的集会次数，那就太好了。
请记住，我的 buffer shapefile 中的每个多边形都有不同的采访日期。
这是我尝试过的方法，但它不起作用：
buffer$count_30 <- GISTools::poly.counts(
rallies[buffer$date - rallies$date > 0 & buffer$date - rallies$date <= 30],
buffer)
我收到以下错误：
Error in `[.data.frame`(x@data, i, j, ..., drop = FALSE) :
undefined columns selected
In addition: Warning messages:
1: In unclass(time1) - unclass(time2) :
longer object length is not a multiple of shorter object length
2: In unclass(time1) - unclass(time2) :
longer object length is not a multiple of shorter object length
完成此任务的正确方法是什么？
我通过使用 sf 包而不是 GISTools 以不同的方式解决了您的问题。该算法简单明了，您可以轻松地将其适应您的 GISTools::poly.counts() 方法：
- 读入 shapefile (
st_read())
- 使用
dplyr 按日期过滤形状文件（确保您有日期对象来创建 windows）
- 找到任何点数据与反弹缓冲区的交集(
st_intersection())
- 获取交集对象的大小(
nrow())
您可能需要调整函数参数以确保它能正确处理真实数据。下面是一个使用模拟数据的示例。
设置并读入数据（注意 stringsAsFactors=F 只是使日期更容易创建；R 版本 4.x 不需要）。
require(tidyverse)
require(magritter) #adds the %<>% operator
require(sf)
require(lubridate)
rally <- st_read(dsn=getwd(),layer='rallies',stringsAsFactors = F)
buff <- st_read(dsn=getwd(),layer='50m_buffer',stringsAsFactors = F)
surv <- st_read(dsn=getwd(),layer='surveys',stringsAsFactors = F)
创建日期对象。
rally %<>% mutate(date=ymd(date))
buff %<>% mutate(date=ymd(date))
surv %<>% mutate(date=ymd(date))
window <- c(ymd('2020-03-27')-30, ymd('2020-03-27')+30)
使用时间过滤数据 window。
buffSub <- buff %>%
filter(date>=window[1] & date<=window[2])
rallySub <- rally %>%
filter(date>=window[1] & date<=window[2])
如果有相交点，则得到数字。
intersectObject <- st_intersection(rallySub, buffSub)
nrow(intersectObject)
或者，如果您想使用自反弹以来的天数或沿着这些线的某些东西，您可以在代表反弹和活动缓冲区之间的时间差的任何点对象中创建新列，并使用这些值进行过滤。
遍历每个集会的日期并获取每个缓冲区的时差。
daysDiff <- data.frame(t(sapply(rally$date, function(d) d-buff$date)))
将这些列添加到数据中并重命名为 buff1、buff2 等
rallyNew <- bind_cols(rally, daysDiff) %>%
rename_with(~gsub('X', 'buff', .x))
使用这些值进行过滤。一次去一列，过滤，并得到与该列关联的缓冲区的交集。
WINDOW=20
for(i in 4:ncol(rallyNew)){
rallySub <- rallyNew %>%
filter(get(unlist(names(rallyNew))[i])<WINDOW &
get(unlist(names(rallyNew))[i])>-WINDOW)
intersectObject <- st_intersection(rallySub, buffSub[i-3,])
print(nrow(intersectObject))
}
另一个答案使用 sf，但这次使用空间连接和 dplyr 进行过滤等
library(tidyverse)
library(sf)
rallies <- read_sf('Downloads/stack_ex_q/rallies.shp')
# Here I don't use the supplied buffer, but make one according to the data
#fifty_buff <- read_sf('Downloads/stack_ex_q/rallies.shp')
surveys <- read_sf('Downloads/stack_ex_q/surveys.shp')
# Transform to a crs using meters as a distance & make date col a proper date
rallies <- st_transform(rallies, crs = 2163) %>%
mutate(date = as.Date(date))
surveys <- st_transform(surveys, crs = 2163) %>%
mutate(date = as.Date(date))
# make a buffer w/ 50 mile radius (80467 meters), not used but useful for visualization
buffer_50mi <- st_buffer(surveys, dist = 80467)
绘制数据以进行快速视觉检查：
library(mapview)
mapview(rallies, col.regions = 'purple') +
mapview(surveys, col.regions = 'black') +
mapview(buffer_50mi, col.regions = 'green')
使用 st_is_within_distance 加入数据，使用 80467m = 50 英里。
joined <- surveys %>%
st_join(rallies, join = st_is_within_distance, 80467)
head(joined)
Simple feature collection with 6 features and 4 fields
geometry type: POINT
dimension: XY
bbox: xmin: 1350401 ymin: -556609 xmax: 1438586 ymax: -455743.1
projected CRS: NAD27 / US National Atlas Equal Area
# A tibble: 6 x 5
id.x date.x geometry id.y date.y
<dbl> <date> <POINT [m]> <dbl> <date>
1 1 2020-04-26 (1350401 -556609) 16 2020-02-19
2 1 2020-04-26 (1350401 -556609) 17 2020-05-12
3 2 2020-03-27 (1438586 -455743.1) 7 2020-02-18
4 2 2020-03-27 (1438586 -455743.1) 15 2020-07-01
5 2 2020-03-27 (1438586 -455743.1) 15 2020-03-28
6 3 2020-06-12 (1352585 -479940.5) 15 2020-07-01
.x 列来自调查 sf 对象，.y 列来自集会 sf 对象。几何形状是从调查 sf 中保留下来的。
使用 dplyr 的过滤器 group_by 和变异，找到您要找的东西。下面的代码以调查点计算 50 英里和 +/- 60 天内的集会作为示例。
joined_60days <- joined %>%
group_by(id.x) %>%
mutate(date_diff = as.numeric(date.x - date.y)) %>%
filter(!is.na(date_diff)) %>% ## remove survey points with no rallies in 50mi/60d
filter(abs(date_diff) <= 60) %>%
group_by(id.x) %>%
count()
head(joined_60days)
Simple feature collection with 4 features and 2 fields
geometry type: POINT
dimension: XY
bbox: xmin: 1268816 ymin: -556609 xmax: 1438586 ymax: -322572.4
projected CRS: NAD27 / US National Atlas Equal Area
# A tibble: 4 x 3
id.x n geometry
<dbl> <int> <POINT [m]>
1 1 1 (1350401 -556609)
2 2 2 (1438586 -455743.1)
3 3 1 (1352585 -479940.5)
4 4 2 (1268816 -322572.4)
快速目视检查：
概览
使用 R，我想根据特定标准（时间window）计算多边形内的点数。
我有以下数据：
- 包含调查访谈日期的地理定位调查数据。因此，我能够准确指出每项调查的时间和地点，并将它们绘制在美国各地。
- 有关美国各地政治集会的地理定位数据。这些还包括日期。
使用 QGIS，我在每个调查受访者周围创建了一组 50 英里的圆形缓冲区。我的目标是计算在采访前的特定时间范围内每个受访者的“缓冲区”内的政治集会次数。 QGIS 中创建的 50 英里缓冲区保留了原始数据的所有变量，包括采访日期。
数据
我使用 QGIS 创建了一些 mock shapefiles 日期和位置以帮助复制。
方法
我正在尝试使用 GISTools::poly.counts 来计算不同时间 windows（30 天、90 天等）内的集会次数。
通常，要计算多边形内的点数，我会简单地使用：
count <- GISTools::poly.counts(rallies, buffer)
这为我提供了每个缓冲区内发生的集会总数，但不允许我指定时间 windows。例如，如果能够统计调查访谈前 30 天和访谈前 90 天缓冲区内的集会次数，那就太好了。
请记住，我的 buffer shapefile 中的每个多边形都有不同的采访日期。
这是我尝试过的方法，但它不起作用：
buffer$count_30 <- GISTools::poly.counts(
rallies[buffer$date - rallies$date > 0 & buffer$date - rallies$date <= 30],
buffer)
我收到以下错误：
Error in `[.data.frame`(x@data, i, j, ..., drop = FALSE) :
undefined columns selected
In addition: Warning messages:
1: In unclass(time1) - unclass(time2) :
longer object length is not a multiple of shorter object length
2: In unclass(time1) - unclass(time2) :
longer object length is not a multiple of shorter object length
完成此任务的正确方法是什么？
我通过使用 sf 包而不是 GISTools 以不同的方式解决了您的问题。该算法简单明了，您可以轻松地将其适应您的 GISTools::poly.counts() 方法：
- 读入 shapefile (
st_read()) - 使用
dplyr按日期过滤形状文件（确保您有日期对象来创建 windows） - 找到任何点数据与反弹缓冲区的交集(
st_intersection()) - 获取交集对象的大小(
nrow())
您可能需要调整函数参数以确保它能正确处理真实数据。下面是一个使用模拟数据的示例。
设置并读入数据（注意 stringsAsFactors=F 只是使日期更容易创建；R 版本 4.x 不需要）。
require(tidyverse)
require(magritter) #adds the %<>% operator
require(sf)
require(lubridate)
rally <- st_read(dsn=getwd(),layer='rallies',stringsAsFactors = F)
buff <- st_read(dsn=getwd(),layer='50m_buffer',stringsAsFactors = F)
surv <- st_read(dsn=getwd(),layer='surveys',stringsAsFactors = F)
创建日期对象。
rally %<>% mutate(date=ymd(date))
buff %<>% mutate(date=ymd(date))
surv %<>% mutate(date=ymd(date))
window <- c(ymd('2020-03-27')-30, ymd('2020-03-27')+30)
使用时间过滤数据 window。
buffSub <- buff %>%
filter(date>=window[1] & date<=window[2])
rallySub <- rally %>%
filter(date>=window[1] & date<=window[2])
如果有相交点，则得到数字。
intersectObject <- st_intersection(rallySub, buffSub)
nrow(intersectObject)
或者，如果您想使用自反弹以来的天数或沿着这些线的某些东西，您可以在代表反弹和活动缓冲区之间的时间差的任何点对象中创建新列，并使用这些值进行过滤。
遍历每个集会的日期并获取每个缓冲区的时差。
daysDiff <- data.frame(t(sapply(rally$date, function(d) d-buff$date)))
将这些列添加到数据中并重命名为 buff1、buff2 等
rallyNew <- bind_cols(rally, daysDiff) %>%
rename_with(~gsub('X', 'buff', .x))
使用这些值进行过滤。一次去一列，过滤，并得到与该列关联的缓冲区的交集。
WINDOW=20
for(i in 4:ncol(rallyNew)){
rallySub <- rallyNew %>%
filter(get(unlist(names(rallyNew))[i])<WINDOW &
get(unlist(names(rallyNew))[i])>-WINDOW)
intersectObject <- st_intersection(rallySub, buffSub[i-3,])
print(nrow(intersectObject))
}
另一个答案使用 sf，但这次使用空间连接和 dplyr 进行过滤等
library(tidyverse)
library(sf)
rallies <- read_sf('Downloads/stack_ex_q/rallies.shp')
# Here I don't use the supplied buffer, but make one according to the data
#fifty_buff <- read_sf('Downloads/stack_ex_q/rallies.shp')
surveys <- read_sf('Downloads/stack_ex_q/surveys.shp')
# Transform to a crs using meters as a distance & make date col a proper date
rallies <- st_transform(rallies, crs = 2163) %>%
mutate(date = as.Date(date))
surveys <- st_transform(surveys, crs = 2163) %>%
mutate(date = as.Date(date))
# make a buffer w/ 50 mile radius (80467 meters), not used but useful for visualization
buffer_50mi <- st_buffer(surveys, dist = 80467)
绘制数据以进行快速视觉检查：
library(mapview)
mapview(rallies, col.regions = 'purple') +
mapview(surveys, col.regions = 'black') +
mapview(buffer_50mi, col.regions = 'green')
使用 st_is_within_distance 加入数据，使用 80467m = 50 英里。
joined <- surveys %>%
st_join(rallies, join = st_is_within_distance, 80467)
head(joined)
Simple feature collection with 6 features and 4 fields
geometry type: POINT
dimension: XY
bbox: xmin: 1350401 ymin: -556609 xmax: 1438586 ymax: -455743.1
projected CRS: NAD27 / US National Atlas Equal Area
# A tibble: 6 x 5
id.x date.x geometry id.y date.y
<dbl> <date> <POINT [m]> <dbl> <date>
1 1 2020-04-26 (1350401 -556609) 16 2020-02-19
2 1 2020-04-26 (1350401 -556609) 17 2020-05-12
3 2 2020-03-27 (1438586 -455743.1) 7 2020-02-18
4 2 2020-03-27 (1438586 -455743.1) 15 2020-07-01
5 2 2020-03-27 (1438586 -455743.1) 15 2020-03-28
6 3 2020-06-12 (1352585 -479940.5) 15 2020-07-01
.x 列来自调查 sf 对象，.y 列来自集会 sf 对象。几何形状是从调查 sf 中保留下来的。
使用 dplyr 的过滤器 group_by 和变异，找到您要找的东西。下面的代码以调查点计算 50 英里和 +/- 60 天内的集会作为示例。
joined_60days <- joined %>%
group_by(id.x) %>%
mutate(date_diff = as.numeric(date.x - date.y)) %>%
filter(!is.na(date_diff)) %>% ## remove survey points with no rallies in 50mi/60d
filter(abs(date_diff) <= 60) %>%
group_by(id.x) %>%
count()
head(joined_60days)
Simple feature collection with 4 features and 2 fields
geometry type: POINT
dimension: XY
bbox: xmin: 1268816 ymin: -556609 xmax: 1438586 ymax: -322572.4
projected CRS: NAD27 / US National Atlas Equal Area
# A tibble: 4 x 3
id.x n geometry
<dbl> <int> <POINT [m]>
1 1 1 (1350401 -556609)
2 2 2 (1438586 -455743.1)
3 3 1 (1352585 -479940.5)
4 4 2 (1268816 -322572.4)
快速目视检查：