[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) <[EMAIL PROTECTED]> Fax-to-email: +44 (0)870 094 0861 Date: 24-Feb-07 Time: 13:45:50 ------------------------------ XFMail ------------------------------ ______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
