Generate random location within specified distance of a given point

Question

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

2. 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 );
}
----------

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

Answer 2

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