[R] random uniform sample of points on an ellipsoid (e.g. WG

[(Ted Harding)]

[Apologies if this is a repeated posting for you. Something seems
 to have gone amiss with my previous attempts to post this reply,
 as seen from my end]

On 22-Feb-07 Roger Bivand wrote:
> On 21 Feb 2007, Russell Senior wrote:
> 
>> 
>> I am interested in making a random sample from a uniform distribution
>> of points over the surface of the earth, using the WGS84 ellipsoid as
>> a model for the earth.  I know how to do this for a sphere, but would
>> like to do better.  I can supply random numbers, want latitude
>> longitude pairs out.
>> 
>> Can anyone point me at a solution?  Thanks very much.
>> 
> 
> http://www.csit.fsu.edu/~burkardt/f_src/random_data/random_data.html
> 
> looks promising, untried.

Hmmm ... That page didn't seem to be directly useful, since
on my understanding of the code (and comments) listed under
"subroutine uniform_on_ellipsoid_map(dim_num, n, a, r, seed, x)"
"UNIFORM_ON_ELLIPSOID_MAP maps uniform points onto an ellipsoid."
in

http://www.csit.fsu.edu/~burkardt/f_src/random_data/random_data.f90

it takes points uniformly distributed on a sphere and then
linearly transforms these onto an ellipsoid. This will not
give unform density over the surface of the ellipsoid: indeed
the example graph they show of points on an ellipse generated
in this way clearly appear to be more dense at the "ends" of
the ellipse, and less dense on its "sides". See:

http://www.csit.fsu.edu/~burkardt/f_src/random_data/
uniform_on_ellipsoid_map.png
[all one line]

Indeed, if I understand their method correctly, in the case
of a horizontal ellipse it is equivalent (modulo rotating
the result) to distributing the points uniformly over a circle,
and then stretching the circle sideways. This will preserve
the vertical distribution (so at the two ends of the major axis
it has the same density as on the circle) but diluting the
horizontal distribution (so that at the two ends of the minor
axis the density isless than on the circle).

I did have a notion about this, but sat on it expecting that
someone would come up with a slick solution -- which hasn't
happened yet.

For the application you have in hand, uniform distribution
over a sphere is a fairly close approximation to uniform
distriobution over the ellipspoid -- but not quite.

But a rejection method, applied to points uniform on the sphere,
can give you points uniform on the ellipsoid and, because of
the close approximation of the sphere to the ellipsoid, you
would not be rejecting many points.

The outline strategy I had in mind (I haven't worked out details)
is based on the following.

Consider a point X0 on the sphere, at radial distance r0 from
the centre of the sphere (same as the centre of the ellipsoid).
Let the radius through that point meet the ellipsoid at a point
X1, at radial distance R1.

Let dS0 be an element of area at X0 on the sphere, which projects
radially onto an element of area dS1 on the ellipsoid. You want
all elements dS1 of equal size to be equally likely to receive
a random point.

Let the angle between the tangent plane to the ellipsoid at X1,
and the tangent plane to the sphere at X0, be phi.

The the ratio of areas dS1/dS0 is R(X0), say, where

  R(X0) = dS1/dS0 = r1^2/(r0^2 * cos(phi))

and the smaller this ratio, the less likely you want a point
u.d. on the sphere to give rise to a point on the ellipsoid.

Now define an acceptance probability P(X0) by

  P(X0) = R(X0)/sup[R(X)]

taking the supremum over X on the sphere. Then sample points X0
unformly on the sphere, accepting each one with probability
P(X0), and continue sampling until you have the number of
points that you need.

Maybe someone has a better idea ... (or code for the above!)

Ted.


--------------------------------------------------------------------
E-Mail: (Ted Harding) <ted.harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 24-Feb-07                                       Time: 13:45:50
------------------------------ XFMail ------------------------------