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

[(Ted Harding)]

Following up to Russell Senior's original query:

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.

As I wrote before:

 "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."

Well, it's "a bit closer than not quite".

I've looked up the WGS84 details, and the ellipsoid has

  major axis a = 6378137.000 metres
  minor axis b = 6356752.314 metres

so b/a = 0.9966471893

As I wrote before:

  "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."


  "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, rejecting each on with probability
  P(X0), and continue sampling until you have the number of
  points that you need."

Without loss of generality, take a = 1 and b as the value of
b'a above, namely

  b = 0.9966471893

Then the minimum value of r1^2/r0^2 is the square of this,
namely 0.99330562 -- this is only 0.67% less than 1.0!

Suppose, therefore, that the ellipsoid has a horizontal
major axis of length a=1, and vertical minor axis of length
b = 0.9966471893, with a circumscribing sphere with the same
centre and radius 1. Consider a planar section containing
the vertical axis. We then have an ellipsoid with horizontal
major axis of length 1 and vertical minor axis of length b.

Let (x,y) be a point of this ellisoid, where

  x = r0*cos(theta), y = r0*sin(theta)

  x^2 + y^2/b^2 = 1

The tangent to the ellipsoid at (x,y) has slope -b^2*cot(theta).

Project (x,y) radially out to the circle. The tangent to the
circle at this point has slope cot(theta).

Hence the angle phi (see above) between the tangent to the
ellipsoid and the tangent to the circle is the difference
between theta and atan(tan(theta)/b^2).

Finally, from the two equations above in (x,y),

  r0 <- 1/sqrt(cos(theta)^2 + sin(theta)^2/b^2)

Therefore we can find out the minimum value of the acceptance
probability above:

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

So we can implement this is some simple R code:

  b <- 0.9966471893
  theta <- (pi/2)*0.001*(0:1000)
  r0 <- 1/sqrt(cos(theta)^2 + sin(theta)^2/b^2)

  phi <- theta - atan(tan(theta)/b^2)
  P <- 1/(r0^2 * cos(phi)) ; P<- P/max(P)

from which min(P) = 0.9933056 -- identical to (b/a)^2 above,
which is not surprising since it occurs where theta = pi/2.

Now, this acceptance probability can be applied to the ellipsoid
as well as to the ellipse, since it involves only the radii and
the angle between the tangent planes; and the the latter is the
same as the angle between the tangents to ellipse and circle.

The question therefore turns on how many points Russell wants
to sample on the ellipsoid of the Earth. If it is a hundred
or two, then there will little perceptible difference between
simply sampling on the sphere and projecting radially onto the
ellipsoid -- maybe 1 or 2 or 3 points should have been rejected,
and quite likely none should have been rejected.

The following code (I hope!) implements the above on the ellipsoid,
testing it for a sample of N=10000, and also prints out the numbers
of points not accepted in each pass.

## Uniform sampling on ellipsoid
b <- 0.9966471893
N <- 10000 ; 
N <- 400
N0 <- N; Accepted <- 0 ;
Lat0 <- numeric(0); Long0 <- numeric(0)
while(Accepted<N0){
  ## Uniform sample on surface of sphere
  x1   <- rnorm(N); x2<-rnorm(N); x3<-rnorm(N)
  requ <- sqrt(x1^2 + x2^2)
  Lat  <- atan2(x3,requ)
  Long <- atan2(x2,x1)
  ## Accept for uniform on ellipsoid
    R0  <- 1/sqrt(cos(Lat)^2 + (sin(Lat)^2)/b^2)
    Phi <- Lat - atan(tan(Lat)/b^2)
    P   <- 1/(R0^2 * abs(cos(Phi))) ; P<- P/max(P)
    Accept   <- (runif(N) <= P)
    Lat0     <- c(Lat0,Lat[Accept])
    Long0    <-c(Long0, Long[Accept])
    Accepted <- Accepted + sum(Accept)
  N <- N0 - Accepted
  print(N)  ## Number not accepted
  ## Another pass if sample size not complete
}

## For instance, on a first pass I got 42 rejected out of 10000,
## then 2 out of the 42, then none out of the final 2.

## Check-plots:
## Polar view:
plot(cos(Long0)*cos(Lat0),sin(Long0)*cos(Lat0),pch="+")
## Equatorial view:
plot(sin(Long0)*cos(Lat0),sin(Lat0),pch="+")
## (These should be similar)

The output of the above code is a set of (Lat0,Long0) points
on the sphere which [should] project radially onto uniformly
distributed points on the surface of the ellipsoid.

I have not gone into the conversion of these (Lat0,Long0)
coordinates on the sphere into latitude and Longitude on
the [WGS80] ellipsoid. This is because there are several
different possible definitions of Latitude (some more
esoteric than others ... ) which can differ by up to 12
minutes of arc. See for instance

  http://en.wikipedia.org/wiki/Latitude

the most "intuitive" ones being "Geographic Latitude" (angle
between equatorial plane and the normal to the reference
ellipsoid) and "Geocentric latitude" (angle between equatorial
plane and radius to point on surface of Earth).

Hoping that helps ...
Ted.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <efh at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 28-Feb-07                                       Time: 23:33:11
------------------------------ XFMail ------------------------------