Обзор

Используя R, я хотел бы подсчитать количество точек внутри многоугольника в соответствии с определенным критерием (временное окно).

У меня есть следующие данные:

1. Географически расположенные данные опроса, которые включают дату опроса. Таким образом, я могу точно определить, когда и где проводилось каждое исследование, и нанести их на карту по всей территории Соединенных Штатов.
2. Географические данные о политических митингах по всей территории Соединенных Штатов. Они также включают дату.

Используя QGIS, я создал набор круговых 50-мильных буферов вокруг каждого респондента опроса. Моя цель — подсчитать количество политических митингов, которые попадают в «буфер» каждого респондента в течение определенного периода времени, предшествующего интервью. 50-километровые буферы, созданные в QGIS, сохраняют все переменные исходных данных, включая дату опроса.

Данные

Используя QGIS, я создал несколько макетных шейп-файлов с датами и местоположениями, чтобы помочь в репликации.

Подход

Я пытаюсь использовать GISTools::poly.counts для подсчета количества митингов в разных временных окнах (30 дней, 90 дней и т.д.).

Как правило, для подсчета количества точек внутри многоугольника я бы просто использовал:

 count <- GISTools::poly.counts(rallies, buffer)
 

Это дает мне общее количество переходов, которые происходят в каждом буфере, но не позволяет мне указывать временные окна. Например, было бы здорово разработать подсчет количества митингов в буфере за 30 дней, предшествующих опросу, а также за 90 дней, предшествующих собеседованию.

Помните, что каждый полигон в моем buffer шейп-файле имеет разную дату собеседования.

Вот что я пробовал, но это не работает:

 buffer$count_30 <- GISTools::poly.counts(
    rallies[buffer$date - rallies$date > 0 amp; 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() методу:

1. Чтение в шейп-файлах ( st_read() )
2. Фильтруйте шейп-файлы по дате с помощью dplyr (убедитесь, что у вас есть объекты даты для создания окон)
3. Найдите пересечение данных любых точек с буфером сбора ( st_intersection() )
4. Получить размер объекта пересечения ( 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)
 

Отфильтруйте данные с помощью временного окна.

 buffSub <- buff %>% 
  filter(date>=window[1] amp; date<=window[2])

rallySub <- rally %>% 
  filter(date>=window[1] amp; date<=window[2])
 

Получите число, если точки пересекаются.

 intersectObject <- st_intersection(rallySub, buffSub)
nrow(intersectObject)
 

Или, если вы хотите использовать дни, прошедшие с момента сбора, или что-то в этом роде, вы можете создать новые столбцы в любом объекте points, которые представляют разницу во времени между сбором и активным буфером, и использовать эти значения для фильтрации.

Перебирайте даты для каждого сбора и получайте разницу во времени для каждого буфера.

 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 amp;
             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 amp; 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 взяты из объекта survey sf, а столбцы .y взяты из объекта rallies sf. Геометрия сохраняется из sf-обзора.

Используя фильтр dplyr, group_by и mutate, найдите то, что вы ищете. В качестве примера приведенный ниже код подсчитывает ралли в пределах 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)
 

Быстрая визуальная проверка: