Random sample of points on a region of the surface of a unit sphere

user_1391441 0
Viewed 430

I need to generate a random sample of points distributed on a region of the surface of a sphere. This answer gives an elegant recipe to sample the entire surface of the sphere:

def sample_spherical(npoints, ndim=3):
    vec = np.random.randn(ndim, npoints)
    vec /= np.linalg.norm(vec, axis=0)
    return vec

(vec returns the (x, y, z) coordinates) which results in:

enter image description here

I need to restrict the sampling to a rectangular region. This region is defined by center values and lengths, for example:

  • center coordinates, eg: (ra_0, dec_0) = (34, 67)
  • side lengths, eg: (ra_l, dec_l) = (2, 1)

where (ra_0, dec_0) are coordinates in the equatorial system, and the lengths are in degrees.

I could use the above function in the following way: call it using a large npoints, transform the returned (x, y, z) values into the equatorial system, reject those outside the defined region, and repeat until the desired number of samples is reached. I feel that there must be a more elegant way of achieving this though.

random geometry python
1 Answers

Assuming that you want to sample points from the 'rectangle' uniformly at random:

  • sample the z-coordinates uniformly at random
  • sample the longitute uniformly at random
import numpy as np

rad = np.pi / 180

def sample_spherical_rect(npoints, ra_0, dec_0, ra_l, dec_l):
    lon_min = (dec_0 - .5 * dec_l) * rad
    lon_delta = dec_l * rad
    z_min = np.sin((ra_0 - .5 * ra_l) * rad)
    z_max = np.sin((ra_0 + .5 * ra_l) * rad)
    z_delta = z_max - z_min
    z = z_min + np.random.rand(npoints) * z_delta
    lat = np.arcsin(z)
    cos_lat = np.cos(lat)
    lon = lon_min + np.random.rand(npoints) * lon_delta
    x = np.cos(lon) * cos_lat
    y = np.sin(lon) * cos_lat
    return np.array([x,y,z])

Explanation: Point of the unit-sphere's surface have the following properties:

  • each of the cartesian corrdinates x,y,z is, on its own, uniform in [-1, 1]

  • given that

    • x = cos(lon) cos(lat)
    • y = sin(lon) cos(lat)
    • z = cos(lat)

    lon is uniform in [0, 2pi].

  • lon and z are statistically independent.

Since a rectangle in the spherical coordinates (lat/lon) is equivalent to a rectangle in the mixed coordinates (z / lon), you just need to sample uniformly at random in (z / lon) from withing your desired ranges and then transform z and lon to x and y.

Further explanation

Since it does not seem to obvious that points (x,y,z) sampled UAR from the unit sphere have a uniform distribution in each of the Cartesian coordinate:

The surface of the slice of the unit sphere bounded by a <= z <= b (with -1 <= a < b <= 1) is given by (can I input latex math stuff in SO?):

A = int_{arcsin(a)}^{arcsin(b)} d[theta] int_0^{2 pi} d[phi] cos(theta) # cos(theta) is the Jacobian
  = 2 pi [sin(theta)]_{theta=arcsin(a)}^{theta=arcsin(b)}
  = 2 pi * (b - a)
user_90962640
Related