6

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.

LoF10
  • 1,907
  • 1
  • 23
  • 64
  • 1
    Measure the distance from all points to the given point (using Pythagoras' theorem), then choose those points with a distance smaller than the desired threshold. – Rodrigo May 19 '19 at 02:48
  • There are various packages that provide circles, but really it can be simply a matter of generating angles (`seq(0,2*pi,len=51)`) and using `sin` and `cos` to turn that into vectors of `x`s and `y`s. You might run into problems, though, as most of this is predicated on having numbers, while you have strings in your data. Another complication: circles of fixed radius are technically more difficult (as are bearing/distance calcs) in lat/lon coords, as pythagorean distances are only approximately correct on a small scale; otherwise, things start to break down (esp at 40deg and more North). – r2evans May 19 '19 at 03:52
  • 1
    @Rodrigo, It's not quite that simple, since lat-lon coordinates are not plane coordinates. You can't just apply the Pythagorean theorem to angles and expect it to all work out. – Cameron Bieganek May 19 '19 at 03:58
  • @clbieganek Yes, you're right, it works as an approximation only. A better option is to calculate the distance in meters, considering the circumference of the Earth at each latitude as the circumference of the Equator multiplied by the cosine of that latitude. Then use the mean latitude for each distance as the basis to measure that distance in meters. That will give a much better approximation. – Rodrigo May 19 '19 at 05:31
  • You should try `destPoint` from `geosphere` package, at it is designed to draw circles on lon/lat coordinates https://stackoverflow.com/a/54638040/7877917 – dieghernan May 30 '19 at 03:47

2 Answers2

7

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.

Cameron Bieganek
  • 7,208
  • 1
  • 23
  • 40
1

The function points_in_circle() from the spatialrisk package handles this.

For example, with your data:

library(spatialrisk)

# Stores 
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)

> spatialrisk::points_in_circle(stores, me$long[1], me$lat[1], radius = 100, lon = long)
# store_id     lat     long distance_m
#        1 40.7505 -73.8456   11.13195
#        2 40.7502 -73.8453   33.70076
mharinga
  • 1,708
  • 10
  • 23