[GIS] Generate random location within specified distance of a given point

circlecoordinatesdistancepointrandom

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:

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

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

Best Answer

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.

Related Solutions

Coordinates Algorithm – Calculating Square Coordinates from Center Point

For this purpose simple approximations are more than good enough. North or south, one degree is about 69 miles but east or west, it is only 69*cos(latitude) miles. Because latitudes do not change much over a ten mile span, you can safely use the cosine of the central latitude of the "square." Therefore the desired coordinates for square vertices at distance r miles from a central location (f,l), given as lat-lon, are computed as

df = r/69        // North-south distance in degrees
dl = df / cos(f) // East-west distance in degrees
{(f-df,l-dl), (f+df,l-dl), (f+df,l+dl), (f-df,l+dl)} // List of vertices

For example, suppose r = 10 miles and the central location is at latitude 50 degrees north, longitude 1 degree west, so that (f,l) = (50,-1) degrees. Then

df = 10/69 = 0.145
dl = 0.145 / cos(50 degrees) = 0.145 / 0.6428 = 0.225
f - df = 50 - 0.145 = 49.855 (southernmost latitude)
f + df = 50 + 0.145 = 50.145 (northernmost latitude)
l - dl = -1 - 0.225 = -1.225 (western longitude)
l + dl = -1 + 0.225 = -0.775 (eastern longitude)

and the coordinates are (49.855,-1.225), (50.145,-1.225), (50.145, -0.775), and (49.855, -0.775) as you march clockwise around the square starting at its southwestern corner.

Don't use this approximation near the poles or for squares larger than a few degrees on a side. Also, depending on limitations of the GIS, some care might be needed around the global cut in longitude, usually taken at +-180 degrees.

Generating Random Locations Nearby – Methods and Tools

This is tricky for two reasons: first, limiting the points to a circle instead of a square; second, accounting for distortions in the distance calculations.

Many GISes include capabilities that automatically and transparently handle both complications. However, the tags here suggest that a GIS-independent description of an algorithm may be desirable.

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

Here is the C code:

Best Answer

Related Solutions

Coordinates Algorithm – Calculating Square Coordinates from Center Point

Generating Random Locations Nearby – Methods and Tools

Related Question