2

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.

Gabriel
  • 40,504
  • 73
  • 230
  • 404

1 Answers1

4

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)
H. Doebler
  • 633
  • 4
  • 9
  • I don't think this is an unbiased estimator. Consider dec_0=45 and dec_l = 90 (that is, the entire northern latitude range). This algorithm would be as likely to produce z=0.9 as it would be to produce z=0.1, despite there being a much wider area around z=0.1. – Sneftel Aug 22 '21 at 14:59
  • @Sneftel No, the area of the unit sphere of the stripe 'z in [0.1, 0.2]' is equal to the area of the stripe 'z in [0.9, 1.0]'. Isnt it? – H. Doebler Aug 22 '21 at 15:06
  • I feel like I must be misunderstanding some aspect of your approach, but as I see it, clearly not. That's why the original approach used a normal distribution, rather than a linear distribution, for each coordinate. – Sneftel Aug 22 '21 at 15:26
  • @Sneftel I added another explanation to my answer. Sampling points (x,y,z) from the whole unit sphere (e.g. using the randn + normalization approad) *results* in random variables x,y,z that are each uniform (aka constant probability density) in [-1, 1], but they are *not independent*. This does not mean that we can simply use three *independent variables* x,y,z to sample points from the unit sphere. But we could sample `phi = 2*pi*rand()` and `theta = arcsin(-1+2*rand())` and return `x,y,z = cos(phi)*cos(theta), sin(phi)*cos(theta), sin(theta)`, effectively sampling `z = -1 + 2*rand()`. – H. Doebler Aug 22 '21 at 15:49
  • That still seems counterintuitive to me, but your math checks out. +1. – Sneftel Aug 22 '21 at 16:07
  • 1
    @Sneftel I totally agree, it is by no means intuitive. I also rubbed my eyes when I first encounetered that property of the sphere. – H. Doebler Aug 22 '21 at 16:57
  • 1
    I think I’m slowly coming around to it. And to answer your earlier question, no… for reasons that still elude me, LaTeX-like markup is supported on many StackExchange sites but not SO. – Sneftel Aug 22 '21 at 17:47