Generate uniform distributed geo points around a geo point (Lat,Long) in a radius

Question

Is there a better way to generate a N numer of Lat/Long coordinates around a Start Point given a distance? This is my approach but it does not creates a uniform distribution.

generate_gps_coordinates_in_circle<-function(start_lat,start_lon,max_distance_in_m,min_distance_in_m){
  
  lat <- (start_lat * pi / 180) #center of the circle
  lon <- (start_lon * pi / 180) #center of the circle
  earth_radius <- 6371000
  distance <- sqrt(runif(1, 0, 1) * (max_distance_in_m ** 2 - min_distance_in_m ** 2) + min_distance_in_m ** 2) #max and min distance from the center
  delta_lat = cos(runif(1) *  pi) * distance / earth_radius
  sign = runif(1, 0, 2) * 2 - 1
  delta_lon = sign * acos(((cos(distance/earth_radius) - cos(delta_lat)) / (cos(lat) * cos(delta_lat + lat))) + 1)
  result <- c(lon = (lon + delta_lon) * 180 / pi, lat =(lat + delta_lat) * 180 / pi)
  return(result)
}

Would it be easier to create a uniform distribution over a square and exclude points greater than a distance `r` from your starting point? — Allan Cameron, Sep 14 '20 at 14:14
I think this is really a probability question rather than a programming question. The difficulty is that circles on a sphere aren't the same as circles in a plane, e.g. if the distance is half the Earth circumference, the length of the boundary is zero; it's just the point on the opposite side of the sphere. — user2554330, Sep 14 '20 at 14:15
Allan this would also work! could you please provide an example? — Andreas, Sep 14 '20 at 14:18
user2554330 hmm could be. I get some how "streched" picture with this solution.. ie points are distributed not in a "round shape" — Andreas, Sep 14 '20 at 14:19
See these posts: CrossValidated [1](https://stats.stackexchange.com/questions/481543/generating-random-points-uniformly-on-a-disk), [2](https://stats.stackexchange.com/questions/79174/how-can-i-generate-uniformly-distributed-points-on-a-circle). — Rui Barradas, Sep 14 '20 at 14:26
Rui i´ve seen those as i researched my problem. But those dont answer my particular case regarding Earth circumferences and geo coordinates as far as i understand. — Andreas, Sep 14 '20 at 14:30
Your question is very unclear. Maybe this: `library(geosphere); p <- geodesic(c(13.408333, 52.518611), azi = runif(100, 0, 360), d = 2000e3); plot(p[, 1:2])` It uses the WGS84 ellipsoid by default (see the documentation). — Roland, Sep 14 '20 at 14:46

score 3 · Accepted Answer · answered Sep 14 '20 at 14:52

You can use the function destPoint from the geosphere package to get the co-ordinates of a given point that is a given bearing and distance from the origin on the Earth's surface. It is simple to get a random bearing by doing runif(1) * 360. However, if you want a uniform density around your point, you will need to take the square root of a uniform random variate as the distance.

We can therefore rewrite your function using destPoint, and also take advantage of vectorization so that we can return a data frame of n points:

generate_gps_coordinates_in_circle <- function(
  start_lat, start_lon, max_distance_in_m, min_distance_in_m, n = 1)
{
  d <- geosphere::destPoint(c(start_lon, start_lat), 360 * runif(n), 
    min_distance_in_m + (max_distance_in_m - min_distance_in_m) * sqrt(runif(n)))
  setNames(as.data.frame(d), c("lon", "lat"))
}

This allows us to do:

df <- generate_gps_coordinates_in_circle(54, 0, 200, 0, n = 1000)

head(df)
#>             lon      lat
#> 1 -0.0010012322 53.99864
#> 2  0.0028912625 53.99959
#> 3  0.0021632655 54.00045
#> 4 -0.0006538438 53.99894
#> 5  0.0014026358 54.00087
#> 6  0.0005882192 54.00172

And we can see the results are satisfactory, producing a circle at small distances:

plot(df$lon, df$lat)

But an appropriately distorted circle at larger distances:

df <- generate_gps_coordinates_in_circle(54, 0, 2000e3, 0, n = 1000)
plot(df$lon, df$lat)

That works well for `max_distance_in_m < 4e6`, but the sphericity of the Earth is pretty noticeable after that. For larger values I'd use `geosphere::randomCoordinates(10*n)`, and reject any points where `geosphere::distGeo()` is bigger than the target distance. — user2554330, Sep 14 '20 at 18:36

Generate uniform distributed geo points around a geo point (Lat,Long) in a radius

1 Answers1