Counting spatial points within a polygon according to a temporal window

Question

Overview

Using R, I would like to count the number of points within a polygon according to a specific criterion (a temporal window).

I have the following data:

Geo-located survey data that include the date of the survey interview. Thus, I am able to pinpoint exactly when and where each survey was conducted and to map them out across the United States.
Geo-located data about political rallies across the United States. These also include the date.

Using QGIS, I created a set of circular 50 mile buffers around each survey respondent. My goal is to count the number of political rallies that falls within each respondent's "buffer" within a specific time frame preceding the interview. The 50 mile buffers created in QGIS retain all variables of the original data, including the date of the interview.

Data

Using QGIS, I created some mock shapefiles with dates and locations to aid in replication.

Approach

I am trying to use GISTools::poly.counts to count the number of rallies within different temporal windows (30 days, 90 days, etc.).

Generally, to count the number of points within a polygon, I would simply use:

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

This gives me the total number of rallies that occur within each buffer, but doesn't allow me to specify temporal windows. For example, it would be great to develop a count of the number of rallies within a buffer for the 30 days preceding the survey interview as well as the 90 days preceding the interview.

Remember, each polygon within my buffer shapefile has a different date of interview.

Here's what I've tried, but it's not working:

buffer$count_30 <- GISTools::poly.counts(
    rallies[buffer$date - rallies$date > 0 & buffer$date - rallies$date <= 30], 
    buffer)

I get the following error:

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

What is the correct way to accomplish this?

user_14781486 · Answer

I approached your problem in a different way by using the sf package instead of GISTools. The algorithm is straightforward and you can easily adapt it to your GISTools::poly.counts() method:

Read in shapefiles (st_read())
Filter shapefiles by date using dplyr (make sure you've got Date objects to create the windows)
Find the intersection of whichever points data with the rally buffer (st_intersection())
Get the size of the intersection object (nrow())

Likely you'll have to tweak the function parameters to make sure it works correctly for the real data. Below is an example using your mock data.

Setup and read in the data (note stringsAsFactors=F just makes the dates easier to create; not necessary for R version 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)

Create the date objects.

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)

Filter the data using the temporal window.

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

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

Get the number if intersecting points.

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

Or if you want to use the days since a rally or something along those lines you can create new columns in whichever points object which represent the time difference between a rally and an active buffer and use those values to filter.

Loop through the dates for each rally and get the time difference with each buffer.

daysDiff <- data.frame(t(sapply(rally$date, function(d) d-buff$date)))

Add those columns to the data and rename with buff1, buff2, etc.

rallyNew <- bind_cols(rally, daysDiff) %>%
  rename_with(~gsub('X', 'buff', .x))

Use those values to filter. Go one column at a time, filter, and get the intersection with buffer associated with that column.

WINDOW=20
for(i in 4:ncol(rallyNew)){
  rallySub <- rallyNew %>% 
    filter(get(unlist(names(rallyNew))[i])-WINDOW)
  intersectObject <- st_intersection(rallySub, buffSub[i-3,])
  print(nrow(intersectObject))
}

user_7547327 · Answer

Another answer using sf, but this time using spatial joins and dplyr for filtering etc.

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)

Plot the data for a quick visual check:

library(mapview)
mapview(rallies, col.regions = 'purple') + 
  mapview(surveys, col.regions = 'black') + 
  mapview(buffer_50mi, col.regions = 'green')

Join the data using st_is_within_distance, using 80467m = 50 miles.

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    
                      
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

The .x columns are from the survey sf object & the .y columns are from the rallies sf object. Geometry is retained from the survey sf.

Using dplyr's filter, group_by, and mutate, find what you're looking for. The code below counts rallies within 50 miles and +/- 60 days by survey point as an example.

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
            
1     1     1   (1350401 -556609)
2     2     2 (1438586 -455743.1)
3     3     1 (1352585 -479940.5)
4     4     2 (1268816 -322572.4)

Quick visual check:

Counting spatial points within a polygon according to a temporal window

Overview

Data

Approach

2 Answers