Filter / subset data between two polygons in R - SF (concentric circle polygons)

Question

Not sure if there's a function i'm missing, but i'm having trouble filtering / checking if point geometries fall between two polygons (concentric circles).

Can you create a mask between two concentric polygons, and then use this to filter out the point geometries that contain the feature data of interest?

I tried subsetting between two polygons using the sf_filter package in R. This did not work.

reproducible code below:

library(sf)
library(tidyverse)
library(sp)

#Create fake data
my.data <- data.frame(replicate(2,sample(-88: -14,100,rep=TRUE))) # Point data
d <- cbind(seq(-180,180,length.out=360),rep(-88,360))
e <- cbind(seq(-180,180,length.out=360),rep(-30,360))

#Project fake data
d = SpatialPoints(cbind(d[,1], d[,2]), proj4string=CRS("+proj=longlat"))
d <- spTransform(d, CRS("+init=epsg:3976"))

e = SpatialPoints(cbind(e[,1], e[,2]), proj4string=CRS("+proj=longlat"))
e <- spTransform(e, CRS("+init=epsg:3976"))

my.data = SpatialPoints(cbind(my.data[,1], my.data[,2]), proj4string=CRS("+proj=longlat"))
my.data <- spTransform(my.data, CRS("+init=epsg:3976"))

d <- sf::st_as_sf(d, coords = c("X1", "X2"),
                  remove = FALSE,
                  crs = st_crs("epsg:3976"))

e <- sf::st_as_sf(e, coords = c("X1", "X2"),
                  remove = FALSE,
                  crs = st_crs("epsg:3976"))

my.data <- sf::st_as_sf(my.data, coords = c("X1", "X2"),
                        remove = FALSE,
                        crs = st_crs("epsg:3976"))
# Create linestrings from circle
d <- d %>% 
  summarise(do_union = FALSE) %>%
  st_cast("LINESTRING")

e <- e %>% 
  summarise(do_union = FALSE) %>%
  st_cast("LINESTRING")

#Join geometries
nst <- rbind(d,e) 

#Create polygon
nst <- nst %>% 
  st_cast("POLYGON")

#Filtering between polygons doesn't return anything
PFz <- st_filter(my.data,nst)

Answer 1

Consider this approach; it builds on three semi random cities in North Carolina (because I love the nc.shp that ships with {sf})

What it does is that it builds two buffers as sf objects, and then constructs two logical vectors - sf::st_contains() for the big circle, and small circle.

Then it is a simple logical operation of checking for points that:

1. are contained within the big circle, and at the same time
2. are not contained within the small circle

Should you want to get more fancy you could run sf::st_difference() on the two buffer objects, and get the mask directly & check for sf::st_contains() only once for the "rim" object.

library(sf)
library(dplyr)

# 3 semi rancom cities in NC (because I *deeply love* the nc.shp file)
cities <- data.frame(name = c("Raleigh", "Greensboro", "Wilmington"),
                     x = c(-78.633333, -79.819444, -77.912222),
                     y = c(35.766667, 36.08, 34.223333)) %>% 
  st_as_sf(coords = c("x", "y"), crs = 4326)

# small buffer - Greensboro will be in; Wilmington not
small_buffer <- cities %>% 
  filter(name == "Raleigh") %>% 
  st_geometry() %>% 
  st_buffer(units::as_units(100, "mile"))

# big buffer - both Greensboro & Wilmington are  in
big_buffer <- cities %>% 
  filter(name == "Raleigh") %>% 
  st_geometry() %>% 
  st_buffer(units::as_units(150, "mile"))

# a visual overview
mapview::mapview(list(big_buffer, small_buffer, cities))

# vector of cities in big buffer
in_big_buffer <- st_contains(big_buffer,
                             cities,
                             sparse = F) %>% 
  t()

# vector of cities in small buffer
in_small_buffer <- st_contains(small_buffer,
                               cities,
                               sparse = F) %>% 
  t()

# cities in concentric circle = in big, and not in small
in_concentric_circle <- in_big_buffer & !in_small_buffer

# check - subset of cities by logical vector
cities %>% 
  filter(in_concentric_circle)

# Simple feature collection with 1 feature and 1 field
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: -77.91222 ymin: 34.22333 xmax: -77.91222 ymax: 34.22333
# Geodetic CRS:  WGS 84
#         name                   geometry
# 1 Wilmington POINT (-77.91222 34.22333)