Find evenly distributed random points on a spherical cap

Question

I have a latitude, longitude and radius of 400m-1000m forming a spherical cap. I need to find a random point on that cap. The points must be evenly distributed over the area.

There is a related question about finding random points in a circle. My first thought was to project the cap on to a Cartesian plane and using the circle algorithm. The radius is small enough that there should be no important level of error.

I'm not sure if projecting and then converting the point back to a lat/lng is the simplest solution or what other possible solutions there are to this problem

Answer 1

You can generate random azimuth in range 0..360 and random distance with sqrt-distribution to provide unifrom distribution

d = maxR * Sqrt(random(0..1))
theta = random(0..1) * 2 * Pi

Then get geopoint coordinates using bearing and distance as described here (Destination point given distance and bearing from start point)

φ2 = asin( sin φ1 ⋅ cos δ + cos φ1 ⋅ sin δ ⋅ cos θ )
λ2 = λ1 + atan2( sin θ ⋅ sin δ ⋅ cos φ1, cos δ − sin φ1 ⋅ sin φ2 )

where   φ is latitude, λ is longitude, θ is the bearing
(clockwise from north), δ is the angular distance d/R; 
d being the distance travelled, R the earth’s radius

Answer 2

As mentioned in the wiki page theta + phi = 90 if phi is a latitude. On the other hand, as r is fixed for all points on the cap, we just need set the value of the theta. Hence, you can pick a random value from 0 to theta value (related to the cap) and define the point by the explained constraints.

Answer 3

For a disc very small compared to the radius of the sphere the lat-long projection will simply be approximately an ellipse except when you are very close to poles.

First compute the stretch at the given latitude:

double k = cos(latitude * PI / 180);

then compute the disc radius in latitude degrees

// A latitude arc-second is 30.87 meters
double R = radius / 30.87 / 3600 * PI / 180;

then compute a uniform random point in a circle

double a = random() * 2 * PI;
double r = R * sqrt(random());

your random point in the disc will be

double random_lat = (latitude*PI/180 + r*cos(a))/PI*180;
double random_longitude = (longitude*PI/180 + (r/k)*sin(a))/PI*180;