# Generate random location within specified distance of a given point

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 *new_long )
{
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 */

/* 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 );
}
----------
``````
• 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. Commented Aug 21, 2013 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. Commented Aug 22, 2013 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? Commented Aug 23, 2013 at 17:55
• @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). Commented Aug 23, 2013 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. Commented Aug 23, 2013 at 19:12

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.

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

For posterity, the formula is:

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
``````