Is there a convenient way to generate points a given distance and direction away from another point on the globe in R?

Question

I'm trying to randomly generate a point within a 10 kilometer radius of a point on the globe in R, and cannot figure out how to in a direct way. The goal is to have an arbitrary point on the globe with coordinates given in WGS84 coordinates, pick a random direction, then a random distance 0-10km and place a new point there. My current solution works, and relies on picking a random distance using runif, then using st_buffer to create a ring of that radius, subtracting the intersection of the ring with a slightly smaller ring to approximate a circle, and then using st_sample to generate a point on the ring. This works, but is incredibly slow, and I need to perform this operation millions of times in my program, so I need to find a more efficient method if possible. The overall goal is to replicate the DHS anonymization algorithm in R on a data set.

I have set up a working means of turning the generated random distances and angles into x and y coordinates in meters that can be added to a point's coordinates to displace it the correct way, but my sample is an sf data frame using WGS84, and the coordinates are therefore not in meters.here are some example coords

I'm not sure how to properly convert the points or translations without having issues involving the curvature of the Earth messing with the distances and angles.

Answer 1

We can use geosphere::destPoint to find the coordinates along some bearing (degrees) and range (meters) from a known point. For an example, the Stack Exchange office is advertised as being in NYC, so we can do

stackexchange <- cbind(-74.00684444589426, 40.708763941777285)
geosphere::destPoint(p=stackexchange, b=145, d=600)
#            lon      lat
# [1,] -74.00277 40.70434

This is vectorized as well. For instance,

• p (the origin coordinates) can be a matrix of 1 or more rows, and if b and d are each length 1 then the output is the same dimensions as the p input, just offset;
• b and d can be length greater than 1 (same lengths), so if p is just 1 row (one coordinate), then it will give you multiple points from that single origin p;
• if all three p, b, and d are the same length (number rows of p, length of b and d), then you get independent bearings/ranges from the points.