Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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:
 *        http://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 );
}
----------
share|improve this question
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
    
Thanks very much for the explanation! –  jetset Aug 23 '13 at 17:36
    
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

3 Answers 3

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.

share|improve this answer

If you are using ArcMap, then consider Point Maker specify your center, number of points and desired distribution. EDIT these sorts of calculations are only useful for projected coordinates not geographic

share|improve this answer
    
This isn't in a GIS system, it's just a C simulation that needs to generate some locations within a radius of a given point. –  jetset Aug 23 '13 at 17:35
    
Randy, please clarify in the post now that you have added the code that this isn't necessarily for a GIS, however, the Python code may provide some guidance –  Dan Patterson Aug 23 '13 at 18:12
    
I don't see any differences between my C and the Python posted in the linked answer. –  jetset Aug 23 '13 at 19:10

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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