I have two tables in different databases I cannot connect directly. Table A contains polygons of regions. Table B contains points with timestamps.

I want to count the points in each polygon for a given timerange. In SQL I would do something like this (simplified reproducible example):

WITH mypolygons(polyid, polygeom) as (
                       6.16527585 46.2437628, 
                       6.16527585 49.62165865, 
                       10.42030488 49.62165865, 
                       10.42030488 46.2437628, 
                       6.16527585 46.2437628
mypoints(pointtime,pointgeom) as (
        '2022-01-01 14:12:00'::timestamp,ST_SetSRID(ST_GeomFromText('Point (7.452982 48.39927093)'),4326)::geometry
timepolygons as (
    from (select 
          generate_series('2022-01-01 14:00:00'::timestamp,'2022-01-01 14:15:00'::timestamp,'5 minute'::interval) as tstamp_start,
          generate_series('2022-01-01 14:05:00'::timestamp,'2022-01-01 14:20:00'::timestamp,'5 minute'::interval) as tstamp_end
     ) as calendar
    cross join 

coalesce(count(mypoints.*),0) as pnt_cnt
left join
mypoints.pointtime between timepolygons.tstamp_start and timepolygons.tstamp_end
group by

However, I cannot use SQL in that case. So I thought of using R with sf package. I already wrote my script which creates the table_polygon, table_points and table_timepolygon as sf dataframes similar to the reproducible example above in SQL.

I know I can count the points in polygons using

(table_polygon$pt_count <- lengths(st_intersects(table_polygon, table_points)))

but how could I add the time-, or more generally spoken, attribute-condition?

The internet is full of questions asking about "count points in polygon", but I could find no examples using an extra attribute or time condition.

As a workaround I thought of using a loop to iterate over each timespan, getting the points falling within and doing the join. But I guess there has to be another, easier way?

1 Answer 1


If you want to go the R way I would suggest this workflow:

  • assign each polygon an unique id / such as the polyid in your example
  • perform sf::st_join(points, polygons) on your points and polygons objects, with - crucially - your points in the first place of the join; it will produce your points object, enriched by the data columns from the polygons object (such as the unique id / polyid).
  • at this step you can drop the geometry (= sf::st_drop_geometry()), and do simple aggregation over your time dimension and polygon id, as both are present in your points object; {dplyr} is your friend here as you are dealing with a regular data frame object

Since an example is worth 1000 of words :) consider this piece of code. It builds on the well known & much loved North Carolina shapefile that ships with {sf}.

It first produces five polygons of dissolved North Carolinas stacked on top of each other, with polygon id's from 1 to 5.

Then it produces three semi-random points, two with time character 1 and one with 3 (or whatever).

The third step is where the action happens: first the points are spatially joined to polygons (it will be a kind of a cross join, since all points are in NC) and then the result is filtered on the secondary criterion - I am using polygon id = time, which makes no sense but is there only for the sake of example. Then it is a question of group by and summarise / functionally equivalent to select polygon id, count(*) from yer data group by polygon id in SQL speak.


# a set of polygons in time
shape <- st_read(system.file("shape/nc.shp", package="sf")) %>%  # included with sf package
  summarise() %>% # single polygon over entire NC
  tidyr::expand_grid(poly_id = 1:5) %>% # sequence of ids from 1 to 5
  st_as_sf(crs = 4267) # this information was lost during the expand_grid (tidyr equivalent of cross join)

# a set of points in time
cities <- data.frame(name = c("Raleigh", "Greensboro", "Wilmington"),
                  x = c(-78.633333, -79.819444, -77.912222),
                  y = c(35.766667, 36.08, 34.223333),
                  time = c(1, 1, 3)) %>%  # time, to be matched with poly_id
  st_as_sf(coords = c("x", "y"), crs = 4326) %>% # CRS of the coordinates
  st_transform(st_crs(shape)) # transorm to CRS of the shape object

# this is the action!   
cities %>% 
  st_join(shape) %>% # spatial join only, without the second criterion (yet)
  st_drop_geometry() %>% # drop geometry (no longer necessary)
  filter(time == poly_id) %>%  # filter over the second condition
  group_by(poly_id) %>%  # group by
  tally(name = "point_totals") # summarise the result
# A tibble: 2 × 2
#  poly_id time_totals
#    <int>       <int>
#1       1           2
#2       3           1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.