library(sf)
nc = st_read(system.file("shape/nc.shp", package="sf"))
p1 = st_sample(nc[1:3, ], 6)
# although coordinates are longitude/latitude, st_intersects assumes that they are planar