Generate random geo coordinates within specific radius from seed point
Asked Answered
C

3

24

I'm using the following function to generate random geo coordinates within a specified radius from a seed point:

function randomGeo(center, radius) {
    var y0 = center.latitude;
    var x0 = center.longitude;
    var rd = radius / 111300;

    var u = Math.random();
    var v = Math.random();

    var w = rd * Math.sqrt(u);
    var t = 2 * Math.PI * v;
    var x = w * Math.cos(t);
    var y = w * Math.sin(t);

    var xp = x / Math.cos(y0);

    return {
        'latitude': y + y0,
        'longitude': xp + x0
    };
}

I do this in a loop, several times, using a 2000m radius and the following seed point:

location: { // Oxford
    latitude: 51.73213,
    longitude: -1.20631
}

I'd expect all of these results to be within 2000m; instead, I'm seeing values upwards of 10000m:

[ { latitude: 51.73256540025445, longitude: -1.3358092771716716 },   // 3838.75070783092
  { latitude: 51.7214165686511, longitude: -1.1644147572878725 },    // 3652.1890457730474
  { latitude: 51.71721400063117, longitude: -1.2082082568884593 },   // 8196.861603477768
  { latitude: 51.73583824510363, longitude: -1.0940424351649711 },   // 5104.820455873758
  { latitude: 51.74017571473442, longitude: -1.3150742602532257 },   // 4112.3279147866215
  { latitude: 51.73496163915278, longitude: -1.0379454413532996 },   // 9920.01459343298
  { latitude: 51.73582333121239, longitude: -1.0939302282840453 },   // 11652.160906253064
  { latitude: 51.72145745285658, longitude: -1.2491630482776055 },   // 7599.550622138115
  { latitude: 51.73036335927129, longitude: -1.3516902043395063 },   // 8348.276271205428
  { latitude: 51.748104753808924, longitude: -1.2669212014250266 },  // 8880.760669882042
  { latitude: 51.72010719621805, longitude: -1.327161328951446 },    // 8182.466715589904
  { latitude: 51.725727610071125, longitude: -1.0691503599266818 } ] // 2026.3687763449955

Given that I (shamelessly!) plagiarized this solution from elsewhere (albeit I've seen several similar implementations), I can't seem to figure out where the math is going wrong.

(Also, in case you want it, this is how I'm calculating the distance. Pretty sure this is correct.)

function distance(lat1, lon1, lat2, lon2) {
    var R = 6371000;
    var a = 0.5 - Math.cos((lat2 - lat1) * Math.PI / 180) / 2 + Math.cos(lat1 * Math.PI / 180) * Math.cos(lat2 * Math.PI / 180) * (1 - Math.cos((lon2 - lon1) * Math.PI / 180)) / 2;
    return R * 2 * Math.asin(Math.sqrt(a));
}
Civilian answered 2/7, 2015 at 18:54 Comment(2)
Did you, by any chance, take it from here? gist.github.com/mkhatib/5641004Prorate
Not there, but that looks the same as mine lolCivilian
R
26

The problem seems to stem from the fact that this is just an inaccurate calculation depending on which center point you are using. Particularly this line:

var xp = x / Math.cos(y0);

Removing this line and changing longitude to

'longitude': x + x0

Seems to keep all of the points within the specified radius, although without this line it seems the points will not completely fill out east to west in some cases.

Anyway, I found someone experiencing a similar issue here with someone elses Matlab code as a possible solution. Depends on how uniformly spread out you need the random points if you wanted to work with a different formula.

Here is a google maps visualization of what's going on with your provided formula:

<!doctype html>
<html>
<head>
    <script type="text/javascript" src="//maps.google.com/maps/api/js?sensor=false"></script>
    <script type="text/javascript" src="//ajax.googleapis.com/ajax/libs/jquery/2.1.4/jquery.min.js"></script>

    <script>
    var distanceLimit = 2000; //in meters
    var numberRandomPoints = 200;
    var mapZoomLevel = 11;
    var locationindex = 0;
    var locations = [
        {'name': 'Oxford, England', 'latitude': 51.73213, 'longitude': -1.20631},
        {'name': 'Quito, Ecuador', 'latitude': -0.2333, 'longitude': -78.5167},
        {'name': 'Ushuaia, Argentina', 'latitude': -54.8000, 'longitude': -68.3000},
        {'name': 'McMurdo Station, Antartica', 'latitude': -77.847281, 'longitude': 166.667942},
        {'name': 'Norilsk, Siberia', 'latitude': 69.3333, 'longitude': 88.2167},
        {'name': 'Greenwich, England', 'latitude': 51.4800, 'longitude': 0.0000},
        {'name': 'Suva, Fiji', 'latitude': -18.1416, 'longitude': 178.4419},
        {'name': 'Tokyo, Japan', 'latitude': 35.6833, 'longitude': 139.6833},
        {'name': 'Mumbai, India', 'latitude': 18.9750, 'longitude': 72.8258},
        {'name': 'New York, USA', 'latitude': 40.7127, 'longitude': -74.0059},
        {'name': 'Moscow, Russia', 'latitude': 55.7500, 'longitude': 37.6167},
        {'name': 'Cape Town, South Africa', 'latitude': -33.9253, 'longitude': 18.4239},
        {'name': 'Cairo, Egypt', 'latitude': 30.0500, 'longitude': 31.2333},
        {'name': 'Sydney, Australia', 'latitude': -33.8650, 'longitude': 151.2094},
    ];
    </script>
</head>
<body>
<div id="topbar">
    <select id="location_switch">
    <script>
        for (i=0; i<locations.length; i++) {
            document.write('<option value="' + i + '">' + locations[i].name + '</option>');
        }
    </script>
    </select>
    <img src="http://google.com/mapfiles/ms/micons/ylw-pushpin.png" style="height:15px;"> = Center
    <img src="https://maps.gstatic.com/mapfiles/ms2/micons/red.png" style="height:15px;"> = No Longitude Adjustment
    <img src="https://maps.gstatic.com/mapfiles/ms2/micons/pink.png" style="height:15px;"> = With Longitude Adjustment (var xp = x / Math.cos(y0);)
</div>

<div id="map_canvas" style="position:absolute; top:30px; left:0px; height:100%; height:calc(100% - 30px); width:100%;overflow:hidden;"></div>

<script>
var markers = [];
var currentcircle;

//Create the default map
var mapcenter = new google.maps.LatLng(locations[locationindex].latitude, locations[locationindex].longitude);
var myOptions = {
    zoom: mapZoomLevel,
    scaleControl: true,
    center: mapcenter
};
var map = new google.maps.Map(document.getElementById('map_canvas'), myOptions);

//Draw default items
var centermarker = addCenterMarker(mapcenter, locations[locationindex].name + '<br>' + locations[locationindex].latitude + ', ' + locations[locationindex].longitude);
var mappoints = generateMapPoints(locations[locationindex], distanceLimit, numberRandomPoints);
drawRadiusCircle(map, centermarker, distanceLimit);
createRandomMapMarkers(map, mappoints);

//Create random lat/long coordinates in a specified radius around a center point
function randomGeo(center, radius) {
    var y0 = center.latitude;
    var x0 = center.longitude;
    var rd = radius / 111300; //about 111300 meters in one degree

    var u = Math.random();
    var v = Math.random();

    var w = rd * Math.sqrt(u);
    var t = 2 * Math.PI * v;
    var x = w * Math.cos(t);
    var y = w * Math.sin(t);

    //Adjust the x-coordinate for the shrinking of the east-west distances
    var xp = x / Math.cos(y0);

    var newlat = y + y0;
    var newlon = x + x0;
    var newlon2 = xp + x0;

    return {
        'latitude': newlat.toFixed(5),
        'longitude': newlon.toFixed(5),
        'longitude2': newlon2.toFixed(5),
        'distance': distance(center.latitude, center.longitude, newlat, newlon).toFixed(2),
        'distance2': distance(center.latitude, center.longitude, newlat, newlon2).toFixed(2),
    };
}

//Calc the distance between 2 coordinates as the crow flies
function distance(lat1, lon1, lat2, lon2) {
    var R = 6371000;
    var a = 0.5 - Math.cos((lat2 - lat1) * Math.PI / 180) / 2 + Math.cos(lat1 * Math.PI / 180) * Math.cos(lat2 * Math.PI / 180) * (1 - Math.cos((lon2 - lon1) * Math.PI / 180)) / 2;
    return R * 2 * Math.asin(Math.sqrt(a));
}

//Generate a number of mappoints
function generateMapPoints(centerpoint, distance, amount) {
    var mappoints = [];
    for (var i=0; i<amount; i++) {
        mappoints.push(randomGeo(centerpoint, distance));
    }
    return mappoints;
}

//Add a unique center marker
function addCenterMarker(centerposition, title) {
    
    var infowindow = new google.maps.InfoWindow({
        content: title
    });
    
    var newmarker = new google.maps.Marker({
        icon: 'http://google.com/mapfiles/ms/micons/ylw-pushpin.png',
        position: mapcenter,
        map: map,
        title: title,
        zIndex: 3
    });
    
    google.maps.event.addListenerOnce(map, 'tilesloaded', function() {
        infowindow.open(map,newmarker);
    });
    
    markers.push(newmarker);
    return newmarker;
}

//Draw a circle on the map
function drawRadiusCircle (map, marker, distance) {
    currentcircle = new google.maps.Circle({
        map: map,
        radius: distance
    });
    currentcircle.bindTo('center', marker, 'position');
}

//Create markers for the randomly generated points
function createRandomMapMarkers(map, mappoints) {
    for (var i = 0; i < mappoints.length; i++) {
        //Map points without the east/west adjustment
        var newmappoint = new google.maps.LatLng(mappoints[i].latitude, mappoints[i].longitude);
        var marker = new google.maps.Marker({
            position:newmappoint,
            map: map,
            title: mappoints[i].latitude + ', ' + mappoints[i].longitude + ' | ' + mappoints[i].distance + 'm',
            zIndex: 2
        });
        markers.push(marker);

        //Map points with the east/west adjustment
        var newmappoint = new google.maps.LatLng(mappoints[i].latitude, mappoints[i].longitude2);
        var marker = new google.maps.Marker({
            icon: 'https://maps.gstatic.com/mapfiles/ms2/micons/pink.png',
            position:newmappoint,
            map: map,
            title: mappoints[i].latitude + ', ' + mappoints[i].longitude2 + ' | ' + mappoints[i].distance2 + 'm',
            zIndex: 1
        });
        markers.push(marker);
    }
}

//Destroy all markers
function clearMarkers() {
    for (var i = 0; i < markers.length; i++) {
        markers[i].setMap(null);
    }
    markers = [];
}

$('#location_switch').change(function() {
    var newlocation = $(this).val();
    
    clearMarkers();

    mapcenter = new google.maps.LatLng(locations[newlocation].latitude, locations[newlocation].longitude);
    map.panTo(mapcenter);
    centermarker = addCenterMarker(mapcenter, locations[newlocation].name + '<br>' + locations[newlocation].latitude + ', ' + locations[newlocation].longitude);
    mappoints = generateMapPoints(locations[newlocation], distanceLimit, numberRandomPoints);

    //Draw default items
    currentcircle.setMap(null);
    drawRadiusCircle(map, centermarker, distanceLimit);
    createRandomMapMarkers(map, mappoints);
});
</script>
</body>
</html>
Ridgeway answered 7/7, 2015 at 22:27 Comment(5)
So far so good, thanks! Curious why it still doesn't fill out east to west? But it does seem to be behaving correctly now.Civilian
My understanding is that your formula would work if the earth was a perfect sphere :) That bit I removed is an attempt to compensate for the earth not being a perfect sphere, but it doesn't work the same every on the globe.Ridgeway
Ahh, you're probably right. Wonder if we could compensate for that? Probably too complicated for what I'm trying to do anyway lolCivilian
I think it gets much more complex. The Haversine Formula seems to be the basis for these calculations. This article on the earth's radius at different latitudes has a clue on how the earth's radius distance number was achieved (it's the mean value of the radius at the equator and the poles). From this I would guess we need a third value along with the coordinates (earth radius at that lat) or could adjust the formula to work correctly for a certain latitude.Ridgeway
Indeed. OK well I'll accept your answer since it solved my immediate problem. if anyone else comes along and wants to contribute a method that supports that 3rd parameter, super! And thanks again ;)Civilian
W
3

You can generate points with a random bearing and distance from the center by moving some distance using vincenty distances (see this stackoverflow answer). In Python, for example, you could use the geopy package.

import random
from geopy import Point
from geopy.distance import geodesic


def generate_point(center: Point, radius: int) -> Point:
    radius_in_kilometers = radius * 1e-3
    random_distance = random.random() * radius_in_kilometers
    random_bearing = random.random() * 360
    return geodesic(kilometers=random_distance).destination(center, random_bearing)


radius = 2000
center = Point(51.73213, -1.20631)
points = [generate_point(center, radius) for _ in range(3000)]

Distances are confirmed with:

assert all(geodesic(center, point).meters <= radius for point in points)
Wartime answered 10/11, 2021 at 17:12 Comment(0)
O
2

Here a simple Vanilla Javascript solution that works like a charm. I want to give credits where it's due and where I found it : https://gist.github.com/fajarlabs/af9e0859fc29b2107bd1797536d2ff2d

/**
* Generates number of random geolocation points given a center and a radius.
* @param  {Object} center A JS object with lat and lng attributes.
* @param  {number} radius Radius in meters.
* @param {number} count Number of points to generate.
* @return {array} Array of Objects with lat and lng attributes.
*/
function generateRandomPoints(center, radius, count) {
  var points = [];
  for (var i=0; i<count; i++) {
    points.push(generateRandomPoint(center, radius));
  }
  return points;
}


/**
* Generates number of random geolocation points given a center and a radius.
* 
* @param  {Object} center A JS object with lat and lng attributes.
* @param  {number} radius Radius in meters.
* @return {Object} The generated random points as JS object with lat and lng attributes.
*/
function generateRandomPoint(center, radius) {
  var x0 = center.lng;
  var y0 = center.lat;
  // Convert Radius from meters to degrees.
  var rd = radius/111300;

  var u = Math.random();
  var v = Math.random();

  var w = rd * Math.sqrt(u);
  var t = 2 * Math.PI * v;
  var x = w * Math.cos(t);
  var y = w * Math.sin(t);

  var xp = x/Math.cos(y0);

  // Resulting point.
  return {'lat': y+y0, 'lng': xp+x0};
}


// Usage Example.
// Generates 100 points that is in a 1km radius from the given lat and lng point.
var randomGeoPoints = generateRandomPoints({'lat':24.23, 'lng':23.12}, 1000, 100);

console.log(randomGeoPoints);
Offshore answered 17/11, 2022 at 16:4 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.