Matlab .
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 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 sample is evenly inside the circle on the azimuthal equidistant projection centered at lat0, lon0. Radial, respectively. azimuthal, the scale for this projection is 1, respectively. g / m Therefore, the range of distortion is r / m, and this is due to the adoption of such points with a probability of t / g.
This code also describes a situation where r0 is about half the circumference of the earth and avoids double sampling, almost antipodal points.
source
share