Logical expressions for spatial geometry and polygon attributes

Question

I want to subset a list of polygons based on (1) whether or not they intersect other polygons in the list in space and (2) an attribute of the polygon that describes time (e.g., when the polygon was created). For example, consider three partially overlapping polygons, a, b, and c created at times t=c(1, 3, 5), respectively. As in the figure below, a and b intersect, b and c intersect, and a and c do not intersect.

I want to subset these polygons to only keep polygons if they intersect another polygon with an earlier time. So for these three polygons, I would exclude a because it intersects b, and b has a later time. I would exclude b because it intersects c, and c has a later time (it intersects a but a has an earlier time; this would not be a reason to exclude it). I would keep c because it intersects b but a has an earlier time.

How can I use both the spatial geometry and the attributes of my polygons to create logical expressions to accomplish this subsetting? I've set up my example problem with the three polygons and time attributes below.

library(sf)

# create polygons
vertices <- rbind(c(1, 1), c(8, 1), c(8, 8), c(1, 8), c(1, 1))
listOfSquares  <- list(a = vertices - 4, b = vertices + 4, c = vertices)
listOfSquares  <-  lapply(listOfSquares, function(x) st_sfc(st_polygon(list(x))) )

# plot polygons
plot(listOfSquares[[1]], xlim = c(-5, 15), ylim = c(-5, 15))
plot(listOfSquares[[2]], add = TRUE)
plot(listOfSquares[[3]], add = TRUE)
text(x=-2.5,y=4.75,"a, t=1")
text(x=-2.5+4,y=4.75+4,"b, t=3")
text(x=-2.5+8,y=4.75+8,"c, t=5")

# add time attribute
listOfSquares[[1]]$time <- 1
listOfSquares[[2]]$time <- 3
listOfSquares[[3]]$time <- 5

Code adapted from this question: R sf: st_intersection on a list of polygons

Answer 1

If you organize the polygons as an sf object you can use st_intersection() to obtain intersections and use dplyr::filter() in combination with dplyr::group_by() to filter down to your desired result.

library(sf)
#> Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2() is TRUE
library(tidyverse)

vertices <- rbind(c(1, 1), c(8, 1), c(8, 8), c(1, 8), c(1, 1))

# The squares b and c are switched in your code. The follwing line corrects that.
listOfSquares  <-
  list(a = vertices - 4, b = vertices, c = vertices + 4)

# Create sf object.
sq_sf <-
  lapply(listOfSquares, function(x)
    st_polygon(list(x))) |>
  st_sfc() |>
  st_sf() |>
  mutate(name = letters[1:3],
         time = c(1, 3, 5))

# Plot with ggplot() to make sure that labels are associated with correct squares.
ggplot(sq_sf) +
  geom_sf() +
  geom_sf_label(aes(label = paste(name, time, sep = ", t=")))

# Get intersections and apply filter rules.
res <-
  st_intersection(sq_sf, sq_sf) |>
  filter(name != name.1) |>
  group_by(name) |>
  filter(all(time.1 < time))
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries

res
#> Simple feature collection with 1 feature and 4 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 5 ymin: 5 xmax: 8 ymax: 8
#> CRS:           NA
#> # A tibble: 1 × 5
#> # Groups:   name [1]
#>   name   time name.1 time.1 st_sfc.lapply.listOfSquares..function.x..st_polygo…¹
#> * <chr> <dbl> <chr>   <dbl>                                            <POLYGON>
#> 1 c         5 b           3                          ((5 5, 5 8, 8 8, 8 5, 5 5))
#> # ℹ abbreviated name:
#> #   ¹​st_sfc.lapply.listOfSquares..function.x..st_polygon.list.x....

# Filter original set to only include the desired squares.
sq_sf |>
  filter(name %in% res$name)
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 5 ymin: 5 xmax: 12 ymax: 12
#> CRS:           NA
#>   st_sfc.lapply.listOfSquares..function.x..st_polygon.list.x.... name time
#> 1                                 POLYGON ((5 5, 12 5, 12 12,...    c    5