PHP Select a random lon/lat within a certain radius

Question

Let's say I have this lon/lat: 33.33333,22.22222

How can I randomly select another lon/lat within an X miles/km radius?

Thanks,

Overlook

It is better to rephrase your original question (#5460561) than post a duplicate. Click the "edit" link below your question to modify your original. – Clergyman

@Clergyman - are you sure? this is a completely different question.. – Overlook

The goal of your question remains the same, this is just a specific part of your original problem. Generally, you will be received better by others if you rephrase your original. It looks like you at least got some traction, but keep that in mind in the future. – Clergyman

I don't have an answer for you but you might find my PHP implementations of geographic calculations might help you github.com/treffynnon/Geographic-Calculations-in-PHP – Carrycarryall

Answer 1

You could use this post to help guide you along:

http://blog.fedecarg.com/2009/02/08/geo-proximity-search-the-haversine-equation/

So with your example, you would just pick a random number between 1 and 10 miles, where 10 is your "within a certain radius".

$longitude = (float) 33.33333;
$latitude = (float) 22.22222;
$radius = rand(1,10); // in miles

$lng_min = $longitude - $radius / abs(cos(deg2rad($latitude)) * 69);
$lng_max = $longitude + $radius / abs(cos(deg2rad($latitude)) * 69);
$lat_min = $latitude - ($radius / 69);
$lat_max = $latitude + ($radius / 69);

echo 'lng (min/max): ' . $lng_min . '/' . $lng_max . PHP_EOL;
echo 'lat (min/max): ' . $lat_min . '/' . $lat_max;

Update:

As Tomalak stated in the comments below, this is working under the assumption that the earth is a sphere rather than a uneven geoid. Because of this, you will get approximations rather than potentially (near)exact results.

Karyogamy

@Tim Nordenfur - are you referring to the 1-10 rand? – Overlook

I wondered what those surprising 69s were doing there, and I think they come from the fact that 1 degree is approximately 69 miles. – Transmigrant

@Spacedman, that is correct... each degree equals approximately 69 miles. – Karyogamy

@MikeLewis: That rather assumes that the Earth is a sphere, when in fact it is an uneven geoid. If you're happy with the approximation then fine, but at least be aware what the margin of error is. In some places it can be quite large. – Gentry

How can I change the code above to support randomization of float values? As my original radius of interest is only 2 miles. – Overlook

@Or W: The probability that the point is 10 miles away is equal to the probability that the point is 1 mile away. If you choose, say, a 1x1 m^2 area close to the original point and one far away, the probability that a random point falls within the closer area is much greater. – Leslie

I've had this function on my system for a long time, I hope this helps you @Or W. codepad.org/bUHx2TWg (Make sure you submit the code to get new results, rather than refreshing the page). – Karyogamy

Thanks for the answer! how would that work for meters? Should the 69 be swapped with 111 045 meters? – Fascicle

Answer 2

@MikeLewis answer is by far a simpler approach, but it only gives you a range of latitude and longitude, and drawing randomly from that might give you points outside the given radius.

The following is a bit more complicated, but should give you 'better' results. (The chances are that isn't necessary, but I wanted to have a go :) ).

As with @MikeLewis' answer the assumption here is that Earth is a sphere. We use this not only in the formulas, but also when we exploit rotational symmetry.

The theory

First we take the obvious approach of picking a random distance $distance (less then $radius miles) and try to find a random point $distance miles away. Such points form a circle on the sphere, and you can quickly convince yourself a straightforward parametrisation of that circle is hard. We instead consider a special case: the north pole.

Points which are a set distance away from the north pole form a circle on the sphere of fixed latitude ( 90-($distance/(pi*3959)*180). This gives us a very easy way of picking a random point on this circle: it will have known latitude and random longitude.

Then we simply rotate the sphere so that our north pole sits at the point we were initially given. The position of our random point after this rotation gives us the desired point.

The code

Note: The Cartesian <--> Spherical co-ordinate transformations used here are different to what is usual in literature. My only motivation for this was to have the z-axis (0,0,1) was pointing North, and the y-axis (0,1,0) pointing towards you and towards the point with latitude and longitude equal to 0. So if you wish to imagine the earth you are looking at the Gulf of Guinea.

/**
 * Given a $centre (latitude, longitude) co-ordinates and a
 * distance $radius (miles), returns a random point (latitude,longtitude)
 * which is within $radius miles of $centre.
 *
 * @param  array $centre Numeric array of floats. First element is 
 *                       latitude, second is longitude.
 * @param  float $radius The radius (in miles).
 * @return array         Numeric array of floats (lat/lng). First 
 *                       element is latitude, second is longitude.
 */
 function generate_random_point( $centre, $radius ){

      $radius_earth = 3959; //miles

      //Pick random distance within $distance;
      $distance = lcg_value()*$radius;

      //Convert degrees to radians.
      $centre_rads = array_map( 'deg2rad', $centre );

      //First suppose our point is the north pole.
      //Find a random point $distance miles away
      $lat_rads = (pi()/2) -  $distance/$radius_earth;
      $lng_rads = lcg_value()*2*pi();


      //($lat_rads,$lng_rads) is a point on the circle which is
      //$distance miles from the north pole. Convert to Cartesian
      $x1 = cos( $lat_rads ) * sin( $lng_rads );
      $y1 = cos( $lat_rads ) * cos( $lng_rads );
      $z1 = sin( $lat_rads );


      //Rotate that sphere so that the north pole is now at $centre.

      //Rotate in x axis by $rot = (pi()/2) - $centre_rads[0];
      $rot = (pi()/2) - $centre_rads[0];
      $x2 = $x1;
      $y2 = $y1 * cos( $rot ) + $z1 * sin( $rot );
      $z2 = -$y1 * sin( $rot ) + $z1 * cos( $rot );

      //Rotate in z axis by $rot = $centre_rads[1]
      $rot = $centre_rads[1];
      $x3 = $x2 * cos( $rot ) + $y2 * sin( $rot );
      $y3 = -$x2 * sin( $rot ) + $y2 * cos( $rot );
      $z3 = $z2;


      //Finally convert this point to polar co-ords
      $lng_rads = atan2( $x3, $y3 );
      $lat_rads = asin( $z3 );

      return array_map( 'rad2deg', array( $lat_rads, $lng_rads ) );
 }

Brahmanism

Answer 3

Pick x1, a number between 0 and x. Pick x2, a number between 0 and x. Your longitude is (1/2)x1 + original longitude and your latitude is (1/2)x2 + original latitude.

Sublittoral

Unfortunately, that won't work either, as an offset of x:x will be outside the radius (distance being x*sqrt(2), or up to 1.4x the distance allowed) – Clergyman

Improvement to my idea: Check the distance of the generated pair from the original coordinates (generated - original) and if it is greater than X, generate coordinates until the distance is smaller than or equal to X. – Sublittoral

Nice thinking. Your edit is better, but now you're either limited to a range of sqrt(2)/2 (only 0.7x the size he was hoping for), or you perform unnecessary duplicate generation until a satisfactory number is created. The only way to properly allow all possibilities is to consider trigonometric functions such as sin/cos. Think of it this way: the area you are creating is a perfect square, while the ideal area is a circle. – Clergyman

I updated it in the comments, not really a good solution since there's probably a faster way than (sort of) brute forcing a solution... But I am still learning kinda basic geometry. – Sublittoral

You are thinking about the problem, and that's the important part. I'm trying to encourage your thinking, not discourage, so don't take it as criticism. :) – Clergyman

Thanks for the help. I really appreciate it. If the radius is only ~0.7x, then it's possible to multiply x by a certain number before doing any other calculations with it, so 0.7x of the x after being multiplied is actually 1x. – Sublittoral

You are correct, but why perform multiple calculations when the most efficient method would be to always generate a valid answer every time? That's the power of using trigonometric functions. It allows you to shift your frame of reference from a cartesian coordinate system (x and y coordinates, or a "grid") to a polar coordinate system (angle and radan coordinates, or a "circle.") Check out en.wikipedia.org/wiki/Polar_coordinate_system to find out more. – Clergyman

Thanks for the link. I don't a lot about trigonometric functions, so I just use what I know. I hope that next year i'll learn more advanced geometry. – Sublittoral

Watched - Very interesting. Thanks again for the link to the video. – Sublittoral

Answer 4

The following Matlab code samples points uniformly on the ellipsoid within a specified distance of a center point.

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.

Tinfoil