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

Alberto Monteiro albmont at centroin.com.br
Fri Mar 2 19:19:50 CET 2007


Russell Senior wrote:
> 
> This is essentially the approach I (and my local helpers) have taken.
> With garden-variety calculus, we have derived a function t(z) that
> maps equal area over the oblate ellipsoid.  We select t from a 
> uniform distribution and then find the corresponding z dimension.  
> The expression is quite hairy, but apparently analytically correct 
> (I haven't found any errors yet), if possibly unstable numerically.
> 
Hmmm... I have just thought about something. The ellipsoid is
"almost" a sphere, so every formula for the sphere should "map"
into a formula for the ellipsoid that can be expressed as a power
series in the eccentricity or the oblateness. Since these parameters
are small, most of these series will converge very fast, and it
might be possible to do most transformations with them.

So, you could have a fast and accurate series that would generate
the random latitudes.

> The derivative with respect to z of area on the ellipsoid at a plane
> of latitude intersecting at a height z above the equator, where the
> ellipsoid has been scaled such that a is the polar radius and the
> equatorial radius is 1, is:
> 
>   dA/dz = 2 * pi * f
> 
> where
> 
>   f = sqrt((z^2 * (1 - a^2))/a^4 + 1)
> 
> You find the function t(z) such that dA/dt is constant.
> 
> Then you select from a uniform distribution and then find the value 
> of z that corresponds.
> 
Yikes! I can't read a formula in ascii notation :-P

Let's parametrize the ellipsoid by the corresponding sphere,
using coordinates lats (latitude of the corresponding sphere)
and lon. So:

x = a cos(lats) cos(lon)
y = a cos(lats) sin(lon)
z = b sin(lats)

But then the surface are between lats and lats + dl will require
a little bit of calculus...

> Russell Senior         ``I have nine fingers; you have ten.''
>
So do Frodo, Sauron and brazilian president Lula :-)

Alberto Monteiro



More information about the R-help mailing list