1

In a simulation (not using any GIS), I need to generate a random location that is within a specified distance of a given location. I found a very helpful post here on this exact topic: How to generate random locations nearby my location?

I see the answer that whuber helpfully provided, but I have two questions on it. First, for ease of reading, here is whuber's answer:


  1. To generate points uniformly, randomly, and independently within a circle of radius r around a location (x0, y0), start by generating two independent uniform random values u and v in the interval [0, 1). (This is what almost every random number generator provides you.) Compute

:

  w = r * sqrt(u)
  t = 2 * Pi * v
  x = w * cos(t) 
  y = w * sin(t)

The desired random point is at location (x+x0, y+y0).

  1. When using geographic (lat,lon) coordinates, then x0 (longitude) and y0 (latitude) will be in degrees but r will most likely be in meters (or feet or miles or some other linear measurement). First, convert the radius r into degrees as if you were located near the equator. Here, there are about 111,300 meters in a degree.

    Second, after generating x and y as in step (1), adjust the x-coordinate for the shrinking of the east-west distances:

    x' = x / cos(y0)
    

    The desired random point is at location (x'+x0, y+y0). This is an approximate procedure. For small radii (less than a few hundred kilometers) that do not extend over either pole of the earth, it will usually be so accurate you cannot detect any error even when generating tens of thousands of random points around each center (x0,y0).


My questions:

First: Can anyone explain what the code is doing?

Specifically:

  w = r * sqrt(u)
  t = 2 * Pi * v
  x = w * cos(t) 
  y = w * sin(t)

Second: When I implemented this, I found it frequently generated locations that slightly exceeded the radius limit. Stepping through the code, everything seemed fine up until the step marked Adjust the x-coordinate for the shrinking of the east-west distances, which is x' = x / cos(y0). That step often caused the longitude offset to get too large, making the resulting longitude jump outside the radius limit. Commenting that out avoids the problem, but I'd really like to understand this better.


Here is the C code:

----------
/*
 * Notes:
 *    Uses information and algorithm from Stackexchange post:
 *        https://gis.stackexchange.com/questions/25877/how-to-generate-random-locations nearby-my-location
 */
int generate_random_location ( float ctr_lat, float ctr_long,
                               float radius_meters, float *new_lat,
                               float *new_long )
{
    float          radius_degrees  = 0; /* radius in degrees */
    float          rand1           = 0; /* a random value between 0 and 1 */
    float          rand2           = 0; /* a random value between 0 and 1 */
    float          w               = 0; /* intermediate value used in formula */
    float          t               = 0; /* intermediate value used in formula */
    float          temp_long_off   = 0; /* intermediate longitude offset (before lat adj) */
    float          lat_off         = 0; /* calculated random latitude offset */
    float          long_off        = 0; /* calculated random longitude offset */
    int            rslt            = 0; /* error indication, if any */


    /* convert input radius from meters to degrees */
    radius_degrees = radius_meters / (float) METERS_PER_DEGREE_AT_EQUATOR;

    /* random() returns a value from 0 to (2**31)-1; divide result by (2**31)-1
       (which is 0xffffffff) to produce a random value between 0 and 1 */
    rand1 = random() / (float) 0xffffffff;
    rand2 = random() / (float) 0xffffffff;

    /* these few lines are Greek to me -- no idea */
    w = radius_degrees * sqrt ( rand1 );
    t = 2 * M_PI * rand2;

    temp_long_off = w * cos ( t );  /* intermediate offset for longitude */
    lat_off       = w * sin ( t );  /* offset for latitude */

    //
    // This adjustment is causing the distance to exceed the radius limit
    //
    /* Adjust longitide offset for the shrinking of east-west distances as
       latitude moves away from the equator */
    long_off = temp_long_off / cos ( ctr_lat );

    /* The random point is at location (lat_off  + center latitude,
                                        long_off + center longitude ) */
    *new_lat  = ctr_lat  + lat_off;
    *new_long = ctr_long + long_off;

    return ( rslt );
}
----------
  • 1
    There are three common reasons for the longitude to appear wrong: (1) you are incorrectly computing the cosine (by mixing up degrees and radians); (2) the radius is so large or the latitude is so close to 90 or -90 degrees that this approximation is no longer any good; (3) there isn't any problem but you are incorrectly visualizing or computing the distances between points. Until you provide more details of what you have done there isn't any way we can reasonably tell which (if any) of these situations you may be in. – whuber Aug 21 '13 at 18:53
  • What the code is doing: "u" is a random distance no greater than radius. Not sure why the SQRT- SQRT of a number from 0 to 1 will still be a number from 0 to 1 but weighted closer to 1. "v" is a random angle, 2 * PI * v converts 0 to 1 to a random angle in radians. The rest of it is just computing position from angle and distance. The x' = x / cos(y0) is a correction if your data is in degrees. It's not needed for projected coordinates. – Russell at ISC Aug 22 '13 at 13:07
  • Let me try to interpret what you said: w = r * sqrt(u) uses the random number u (0-1) and the radius in degrees to create a random distance within the radius. The SQRT biases the distance to be closer to the maximum. That might be a problem in my simulation. 2. t = 2 * Pi * v uses the random number v (0-1) to create a random angle in radians 3. Since t is then in degrees, the correction x' = x / cos(y0) isn't needed. 3. Did I get that right? – jetset Aug 23 '13 at 17:55
  • 1
    @Russell The sqrt unbiases the distance: without it, too many points would appear too close to the origin. This is a consequence of the Pizza principle, which implies that fully half the points would appear within just one-quarter of the area if you did not take the square root. Randy, the cosine correction is needed for geographic (lat, lon) coordinates. It is not needed for projected coordinates (but then the projection needs to be sufficiently accurate). – whuber Aug 23 '13 at 18:26
  • whuber, thanks for the clarifications. Can you shed any light on why the cosine correction shifts the longitude such that it often slightly exceeds the input radius? Even when I tried as a hack to reduce the radius by 25% beforehand. – jetset Aug 23 '13 at 19:12
1

If you are serious about sampling uniformly within a distance of a particular point then you need account for the curvature of the ellipsoid. This is pretty easy to do using a rejection technique. The following Matlab code implements this.

function [lat, lon] = geosample(lat0, lon0, r0, n)
% [lat, lon] = geosample(lat0, lon0, r0, n)
%
% Return n points on the WGS84 ellipsoid within a distance r0 of
% (lat0,lon0) and uniformly distributed on the surface.  The returned
% lat and lon are n x 1 vectors.
%
% Requires Matlab package
%  http://www.mathworks.com/matlabcentral/fileexchange/39108

  todo = true(n,1); lat = zeros(n,1); lon = lat;

  while any(todo)
    n1 = sum(todo);
    r = r0 * max(rand(n1,2), [], 2);  % r = r0*sqrt(U) using cheap sqrt
    azi = 180 * (2 * rand(n1,1) - 1); % sample azi uniformly
    [lat(todo), lon(todo), ~, ~, m, ~, ~, sig] = ...
        geodreckon(lat0, lon0, r, azi);
    % Only count points with sig <= 180 (otherwise it's not a shortest
    % path).  Also because of the curvature of the ellipsoid, large r
    % are sampled too frequently, by a factor r/m.  This following
    % accounts for this...
    todo(todo) = ~(sig <= 180 & r .* rand(n1,1) <= m);
  end
end

This code samples uniformly within a circle on the azimuthal equidistant projection centered at lat0, lon0. The radial, resp. azimuthal, scale for this projection is 1, resp. r/m. Hence the areal distortion is r/m and this is accounted for by accepting such points with a probability m/r.

This code also accounts for the situation where r0 is about half the circumference of the earth and avoids double sampling nearly antipodal points.

0

You can use this formula, with a fixed distance d and random value of tc: http://williams.best.vwh.net/avform.htm#LL

All parameters are radians.

edit: This link is dead. Here is an archived version: https://web.archive.org/web/20161209044600/http://williams.best.vwh.net/avform.htm#LL

For posterity, the formula is:

Lat/lon given radial and distance

A point {lat,lon} is a distance d out on the tc radial from point 1 if:

 lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
 IF (cos(lat)=0)
    lon=lon1      // endpoint a pole
 ELSE
    lon=mod(lon1-asin(sin(tc)*sin(d)/cos(lat))+pi,2*pi)-pi
 ENDIF

This algorithm is limited to distances such that dlon extend around less than one quarter of the circumference of the earth in longitude. A completely general, but more complicated algorithm is necessary if greater distances are allowed:

 lat =asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
 dlon=atan2(sin(tc)*sin(d)*cos(lat1),cos(d)-sin(lat1)*sin(lat))
 lon=mod( lon1-dlon +pi,2*pi )-pi

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.