How to generate random geocoordinates using random module

The problem when using uniform distributions for both the latitude and the longitude is that physically, the latitude is NOT uniformly distributed.

So if you plan to use these random points for something like some statiscal averaging computation, or a physics Monte-Carlo simulation, the results will risk being incorrect.

And if you plot a graphical representation of the “uniform” random points, they will seem to cluster in the polar regions.

To picture that, consider on planet Earth the zone that lies between 89 and 90 degrees of latitude (North). The length of a degree of latitude is 10,000/90 = 111 km. That zone is a circle of radius 111 km, centered around the North Pole. Its area is about 3.14 * 111 * 111 ≈ 39,000 km²

On the other hand, consider the zone that lies between 0 and 1 degree of latitude. This is a strip whose length is 40,000 km (the Equator) and whose width is 111 km, so its area is 4.44 millions km². Much larger than the polar zone.

A simple algorithm:

A possibility is to use Gaussian-distributed random variables, as provided by the Python library. If we build a 3D vector whose 3 components have Gaussian distributions, the overall probability distribution is like exp(-x²) * exp(-y²) * exp(-z²) but this is the same thing as exp(-(x² + y² + z²)) or exp(-r²), where r is the distance from the origin.

So these vectors have no privileged direction. Once normalized to unit length, they are uniformly distributed on the unit sphere. They solve our problem with the latitude distribution.

The idea is implemented by the following Python code:

import  math
import  random

def randlatlon1():
    pi = math.pi
    cf = 180.0 / pi  # radians to degrees Correction Factor

    # get a random Gaussian 3D vector:
    gx = random.gauss(0.0, 1.0)
    gy = random.gauss(0.0, 1.0)
    gz = random.gauss(0.0, 1.0)

    # normalize to an equidistributed (x,y,z) point on the unit sphere:
    norm2 = gx*gx + gy*gy + gz*gz
    norm1 = 1.0 / math.sqrt(norm2)
    x = gx * norm1
    y = gy * norm1
    z = gz * norm1

    radLat = math.asin(z)      # latitude  in radians
    radLon = math.atan2(y,x)   # longitude in radians

    return (round(cf*radLat, 5), round(cf*radLon, 5))

A sanity check:

Euclidean geometry provides a formula for the probability of a spherical zone defined by minimal/maximal latitude and longitude. The corresponding Python code is like this:

def computeProbaG(minLat, maxLat, minLon, maxLon):
    pi  = math.pi
    rcf = pi / 180.0  # degrees to radians Correction Factor
    lonProba = (maxLon - minLon) / 360.0
    minLatR  = rcf * minLat
    maxLatR  = rcf * maxLat
    latProba = (1.0/2.0) * (math.sin(maxLatR) - math.sin(minLatR))
    return (lonProba * latProba)

And we can also compute an approximation of that same probability by random sampling, using the random points provided by a function such as randlatlon1, and counting what percentage of them happen to fall within the selected zone:

def computeProbaR(randlatlon, ranCount, minLat, maxLat, minLon, maxLon):
    norm  = 1.0 / ranCount
    pairs = [randlatlon()  for i in range(ranCount)]
    acceptor = lambda p: ( (p[0] > minLat)  and  (p[0] < maxLat)     and
                           (p[1] > minLon)  and   (p[1] < maxLon) )
    selCount = sum(1 for p in filter(acceptor, pairs))
    return (norm * selCount)

Equipped with these two functions, we can check for various geometric parameter sets that the geometric and probabilistic results are in good agreement, with ranCount set to one million random points:

ranCount = 1000*1000
print (" ")
probaG1 = computeProbaG(                       30, 60, 45, 90)
probaR1 = computeProbaR(randlatlon1, ranCount, 30, 60, 45, 90)
print ("probaG1 = %f" % probaG1)
print ("probaR1 = %f" % probaR1)
print (" ")
probaG2 = computeProbaG(                        10, 55, -40, 160)
probaR2 = computeProbaR(randlatlon1, ranCount,  10, 55, -40, 160)
print ("probaG2 = %f" % probaG2)
print ("probaR2 = %f" % probaR2)
print (" ")

Execution output:

$ python3 georandom.py

probaG1 = 0.022877
probaR1 = 0.022852
 
probaG2 = 0.179307
probaR2 = 0.179644

$

So the two sort of numbers appears to agree reasonably here.

Addendum:

For the sake of completeness, we can add a second algorithm which is less intuitive but derives from a wider statistical principle.

To solve the problem of the latitude distribution, we can use the Inverse Transform Sampling theorem. In order to do so, we need some formula for the probability of the latitude to be less than an arbitrary prescribed value, φ.

The region of the unit 3D sphere whose latitude is less than a given φ is known as a spherical cap. Its area can be obtained thru elementary calculus, as described here for example.

The spherical cap area is given by formula: A = 2π * (1 + sin(φ)) The corresponding probability can be obtained by dividing this area by the overall area of the unit 3D sphere, that is 4π, corresponding to φ = φ_max = π/2. Hence:

p = Proba{latitude < φ} = (1/2) * (1 + sin(φ))

Or, conversely:

φ = arcsin (2*p - 1)

From the Inverse Transform Sampling theorem, a fair sampling of the latitude (in radians) is obtained by replacing the probability p by a random variable uniformly distributed between 0 and 1. In Python, this gives:

lat = math.asin(2*random.uniform(0.0, 1.0) - 1.0)

As for the longitude, this is an independent random variable that is still uniformly distributed between -π and +π (in radians). So the overall Python sampler code is:

def randlatlon2r():
    pi = math.pi
    cf = 180.0 / pi  # radians to degrees Correction Factor
    u0 = random.uniform(0.0, 1.0)
    u1 = random.uniform(0.0, 1.0)
    radLat = math.asin(2*u0 - 1.0)  # angle with Equator   - from +pi/2 to -pi/2
    radLon = (2*u1 - 1) * pi        # longitude in radians - from -pi to +pi 
    return (round(radLat*cf,5), round(radLon*cf,5))

This code has been found to pass successfully the sanity check as described above.