simple approach to determine if points intersect overlapping polygons

I have a simple dataframe with latitude and longitude coordinates (I've provided a subset below). I'd like to determine if any of the points fall within one or both polygons (the second polygon completely overlaps the first).

# subsample of points
pts <- structure(list(Longitude = c(-26.775, -27.3167, -30.2717, -30.7267, -37.8683, -38.2733, 
                                    -44.3033, -44.8233, -57.42, -57.8717, -52.235, -52.54, 
                                    -54.2317, -54.6883, -63.6033, -65.415, -66.31, -66.7533), 
                      Latitude = c(62.5183, 62.3, 60.8367, 60.59, 56.6383, 56.3917, 52.8, 52.45, 
                                   45.56, 45.46, 47.3917, 47.1317, 46.26, 46.16, 44.29, 43.2633, 
                                   43.0983, 43.0267)), 
                 row.names = c(NA, -18L), class = c("tbl_df", "tbl", "data.frame"))

# points as a spatial object
sp <- SpatialPoints(pts, proj4string = CRS("+proj=longlat +datum=WGS84"))

# polygon A & B
pA <- st_polygon(list(as.matrix(data.frame(lon = c(-40, -40, -64, -64, -54, -54, -40),
                                              lat = c(45, 61, 61, 51, 51, 45, 45))))) %>%
  st_sfc(crs = st_crs("+proj=longlat +datum=WGS84"))
pB<- st_polygon(list(as.matrix(data.frame(lon = c(-6, -6, -50, -50, -64, -64, -54, -54, -6),
                                                 lat = c(45, 66, 66, 70, 70, 51, 51, 45, 45))))) %>%
  st_sfc(crs = st_crs("+proj=longlat +datum=WGS84"))

I know I can use st_intersects to determine if points overlap in a single polygon, but I don't know how to expand this if when you have multiple overlapping polygons. Ideally, the end result would be a variable that returns A if the point is just in A, B if the point is part of the larger B polygon that does not overlap A, and NA if the point is not in either polygon.

I find it much easier to follow code that uses tidyverse (rather than datatable), so I'd prefer a solution with those functions.

Solution

If you have collected your polygons into sf object and find point intersects with that, st_intersects() will return a "sparse geometry binary predicate list", it includes a list of intersecting geometry indices for each point:

Sparse geometry binary predicate list of length 18, where the predicate was `intersects'
first 10 elements:
 1: 1
 2: 1
 3: 1
 4: 1
 5: 1
 6: 1
 7: 1, 2
 8: 1, 2
 9: (empty)
 10: (empty)

If your polygons are ordered correctly in sf, last item for every point will be the polygon with highest row number. Well, if you just want know if point intersects with 1, 2, (,n) or no polygons, all you really need is to get the lengths() of that list.

library(sf)
library(sp)
library(dplyr)
library(ggplot2)

### move to sf objects, assume that there are no partial overlaps,
### arrange polygons by their size
ab_sf <- st_sf(geometry = c(pA, pB), row.names = c("pA", "pB")) %>% 
  tibble::rownames_to_column("poly") %>% 
  mutate(size_rank = rank(desc(st_area(geometry)))) %>%  
  arrange(size_rank)
ab_sf
#> Simple feature collection with 2 features and 2 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -64 ymin: 45 xmax: -6 ymax: 70
#> Geodetic CRS:  +proj=longlat +datum=WGS84
#>   poly size_rank                       geometry
#> 1   pB         1 POLYGON ((-6 45, -6 66, -50...
#> 2   pA         2 POLYGON ((-40 45, -40 61, -...

points_sf <- st_as_sf(sp_)

### find intersects, returned list will include ordered list of intersecting indices
### i.e (1) for only pA; (1,2) for pA and pB; (empty) for no intersect,
### extract last item from each to get index of smallest intersecting polygon
isect_last_idx <- st_intersects(points_sf, ab_sf) %>% sapply(last)
isect_last_idx
#>  [1]  1  1  1  1  1  1  2  2 NA NA  2  2 NA NA NA NA NA NA

## map indeces to polygon names
points_sf$poly_name <- ab_sf[["poly"]][isect_last_idx]
points_sf
#> Simple feature collection with 18 features and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -66.7533 ymin: 43.0267 xmax: -26.775 ymax: 62.5183
#> Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs
#> First 10 features:
#>                    geometry poly_name
#> 1   POINT (-26.775 62.5183)        pB
#> 2     POINT (-27.3167 62.3)        pB
#> 3  POINT (-30.2717 60.8367)        pB
#> 4    POINT (-30.7267 60.59)        pB
#> 5  POINT (-37.8683 56.6383)        pB
#> 6  POINT (-38.2733 56.3917)        pB
#> 7     POINT (-44.3033 52.8)        pA
#> 8    POINT (-44.8233 52.45)        pA
#> 9      POINT (-57.42 45.56)      <NA>
#> 10   POINT (-57.8717 45.46)      <NA>

ab_sf %>% 
  arrange(size_rank) %>% 
  ggplot() +
  geom_sf(aes(fill = poly, linetype = poly, linewidth = poly), alpha = .8) +
  scale_linewidth_manual(values = c(0,2)) +
  geom_sf(data = points_sf,  aes(shape = poly_name), size =  5) +
  scale_shape_manual(na.value = 1, values = c(9,13)) +
  theme_minimal()

Used data:


pts <- structure(list(Longitude = c(-26.775, -27.3167, -30.2717, -30.7267, -37.8683, -38.2733, 
                                    -44.3033, -44.8233, -57.42, -57.8717, -52.235, -52.54, 
                                    -54.2317, -54.6883, -63.6033, -65.415, -66.31, -66.7533), 
                      Latitude = c(62.5183, 62.3, 60.8367, 60.59, 56.6383, 56.3917, 52.8, 52.45, 
                                   45.56, 45.46, 47.3917, 47.1317, 46.26, 46.16, 44.29, 43.2633, 
                                   43.0983, 43.0267)), 
                 row.names = c(NA, -18L), class = c("tbl_df", "tbl", "data.frame"))

# points as a spatial object
sp_ <- SpatialPoints(pts, proj4string = CRS("+proj=longlat +datum=WGS84"))

# polygon A & B
pA <- st_polygon(list(as.matrix(data.frame(lon = c(-40, -40, -64, -64, -54, -54, -40),
                                           lat = c(45, 61, 61, 51, 51, 45, 45))))) %>%
  st_sfc(crs = st_crs("+proj=longlat +datum=WGS84"))
pB<- st_polygon(list(as.matrix(data.frame(lon = c(-6, -6, -50, -50, -64, -64, -54, -54, -6),
                                          lat = c(45, 66, 66, 70, 70, 51, 51, 45, 45))))) %>%
  st_sfc(crs = st_crs("+proj=longlat +datum=WGS84"))

^{Created on 2023-05-25 with reprex v2.0.2}