How to randomly select a point on the surface of the Earth?

Question

As the title says, how to randomly select a coordinate on the surface of the Earth without bias? The simple solution of doing rand(-90,90),rand(-180,180) will favor polar regions over equatorial regions.

Answer 1

It sounds like you're looking for a uniform random distribution on a sphere.

In spherical coordinates:

θ = 2π * rand(0,1)
Φ = arccos(1 - 2*rand(0,1))

Answer 2

Archimedes proved a nifty theorem that the area of a region on a sphere is equal to the area of the horizontal projection of that region onto the curved face of a circumscribing cylinder (and hence in particular, the area of the whole sphere is equal to the area of the whole curved face of the circumscribed cylinder). You can convince yourself of this by drawing little "rectangles" (bounded by lines of latitude and longitude) and proving it for those, then observing (informally) that at the limit everything is approximated by small enough rectangles, or more formally that the area of a region is defined to be some integral, i.e. (thanks to Lebesgue) a limit of a series of sums of rectangles.

So, if you choose a random uniformly-distributed point on the curved face of the cylinder (easy since it's isometric to a rectangle) and then project that point horizontally back onto the sphere, your distribution will have the property that the probability of a point being in any region is proportional to the area of the region.

This means choose the longitude uniformly in the range (-pi, pi], and set the latitude to arccos(y) - pi/2, where y is uniformly distributed in the range (-1,1].

You'll still have a small bias thanks to the (lack of) precision of whatever floating-point type you're using for arccos. I'm not sure whether that will be noticeable, but if so and if you want to fix that you can make your y-values more precise close to +/-1 (in effect by generating more bits after the decimal point). Or rather, split the range into two, generate extra bits close to 0, then use an accurate arccos_1_minus_x function.