R - How do I draw a radius around a point and use that result to filter other points?

Question

I am looking to draw a radius around a lat long point, then use that buffer to filter other points that fit within it. For example:

#stores datasets
stores = data.frame(store_id = 1:3,
                    lat = c("40.7505","40.7502","40.6045"),
                    long = c("-73.8456","-73.8453","-73.8012")
                    )

#my location
me  = data.frame(lat = "40.7504", long = "-73.8456")

#draw a 100 meter radius around me 


#use the above result to check which points in dataset stores are within that buffer

Not sure how to approach this. I've worked with over before to intersect points and polygons but not sure how to run a similar scenario on lone points.

Answer 1

You could hypothetically try to do geometric calculations on the surface of a sphere or ellipsoid, but usually when one is performing geometric map operations, one uses a map projection which projects the lon-lat coordinates onto a flat surface.

Here is how this can be done using the sf package. First, create your points in lon-lat coordinates:

library(sf)

lat <- c(40.7505, 40.7502, 40.6045)
lon <- c(-73.8456, -73.8453, -73.8012)

stores <- st_sfc(st_multipoint(cbind(lon, lat)), crs = 4326)

me <- st_sfc(st_point(c(-73.8456, 40.7504)), crs = 4326)

The crs = 4326 argument specifies the EPSG code for the lon-lat coordinate system. Next we need to pick a map projecton. For this example I will use UTM zone 18, which contains the above points:

stores_utm <- st_transform(stores, "+proj=utm +zone=18")
me_utm     <- st_transform(me, "+proj=utm +zone=18")

Now we can buffer the point representing ourself by 100 meters to produce a circle with radius 100 meters:

circle <- st_buffer(me_utm, 100)

Now, we're almost ready to use a geometric predicate to test which points are in the circle. However, stores_utm is currently a MULTIPOINT, so geometric predicates will treat it like one geometric entity. We can fix this by casting stores_utm as a POINT, which will give us a collection of three separate points:

stores_utm_column <- st_cast(stores_utm, "POINT")
stores_utm_column
# Geometry set for 3 features 
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 597453 ymin: 4495545 xmax: 601422.3 ymax: 4511702
# epsg (SRID):    32618
# proj4string:    +proj=utm +zone=18 +ellps=WGS84 +units=m +no_defs
# POINT (597453 4511702)
# POINT (597478.7 4511669)
# POINT (601422.3 4495545)

Now we can test which points are in the circle:

> st_contains(circle, stores_utm_column, sparse = FALSE)
#      [,1] [,2]  [,3]
# [1,] TRUE TRUE FALSE

which shows that the first two points are in the circle and the third one is not.

Of course, every map projection introduces some distortions. Your choice of projection will depend on the nature of your problem.