I would like to generate coordinates within a given radius of given coordinates. I wrote a little script, that most of the time generate valid values:

const coordinate = [Math.random() * 180 - 90, Math.random() * 360 - 180]; // [lat, long]
const angleRadians = Math.PI / 4; // "random" angle in ° radians
const totalDistance = 100; // max distance in km

// First calculate the new latitude
const kmPerDegreeLatitude = 111; // in km/°
const distanceLatitude = totalDistance / kmPerDegreeLatitude; // in °

const offsetLatitude = Math.sin(angleRadians) * distanceLatitude;
const newLatitude = coordinate[0]   offsetLatitude;

// This looks wrong to me, but reduces the likelihood of errors
const remainingDistance = Math.sqrt(
  Math.pow(totalDistance, 2) -
    Math.pow(offsetLatitude * kmPerDegreeLatitude, 2)
);

// Now calculate the longitude
const kmPerDegreeLongitude =
  Math.abs(Math.cos(degreesToRadians(newLatitude))) * 111; // in km/°
const distanceLongitude = remainingDistance / kmPerDegreeLongitude; // in °

const offsetLongitude = Math.cos(angleRadians) * distanceLongitude;
const newLongitude = coordinate[1]   offsetLongitude;

// Done
const newCoordinate: [latitude: number, longitude: number] = [
  newLatitude,
  newLongitude,
];

But when I use the following code to check the distance I sometimes end up above the allowed distance:

function haversine(
  latitude1: number,
  longitude1: number,
  latitude2: number,
  longitude2: number,
) {
  const EQUATORIAL_EARTH_RADIUS = 6378.137;
  const distanceLatitude = degreesToRadians(latitude2 - latitude1);
  const distanceLongitude = degreesToRadians(longitude2 - longitude1);
  const a =
    Math.sin(distanceLatitude / 2) * Math.sin(distanceLatitude / 2)  
    Math.cos(degreesToRadians(latitude1)) *
      Math.cos(degreesToRadians(latitude2)) *
      Math.sin(distanceLongitude / 2) *
      Math.sin(distanceLongitude / 2);
  const distance =
    EQUATORIAL_EARTH_RADIUS * 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));

  return distance;
}

Any idea, where I messed up? I assume it's in the coordinate generation, but it might be in the verification step. Or are there maybe way easier ways to generate the nearby coordinates?

Please note, that I'm not able to use existing libraries for that, because the generated coordinates should be reproducible, but I omitted the related code for simplicity.

CodePudding user response：

Well, I tried to overcomplicate things and the float imprecision in JS stacked up sometimes to 20% due to stacked calls.

I simplified the calculation now to assume we live on a perfect sphere:

const angleRadians = randomAngle(); // in ° radians

const errorCorrection = 0.995; // avoid float issues
const distanceInKm = randomDinstance(100) * errorCorrection; // in km

const distanceInDegree = distanceInKm / kmPerDegree; // in °

const newCoordinate: [latitude: number, longitude: number] = [
  coordinate[0]   Math.sin(angleRadians) * distanceInDegree,
  coordinate[1]   Math.cos(angleRadians) * distanceInDegree,
];

// Box latitude [-90°, 90°]
newCoordinate[0] = newCoordinate[0] % 180;
if (newCoordinate[0] < -90 || newCoordinate[0] > 90) {
  newCoordinate[0] = Math.sign(newCoordinate[0]) * 180 - newCoordinate[0];
  newCoordinate[1]  = 180;
}
// Box longitude [-180°, 180°]
newCoordinate[1] = (((newCoordinate[1] % 360)   540) % 360) - 180;

CodePudding user response：

Assuming that you can locally approximate the surface of the earth as a plane, then the following function gives you the latitude and longitude of a point P1, given the latitude and longitude of a point P0, the linear distance |P1-P0| and the angle formed by P1-P0 with the local parallel. The local earth radius R is also necessary as input.

function incrementCoordinates (long0, lat0, dist, angle, R) {

    // Calculate the distance component along the parallel
    dist_x = dist * Math.cos(angle / 180 * Math.PI)

    // Calculate the distance component along the meridian
    dist_y = dist * Math.sin(angle / 180 * Math.PI)

    // Calculate the new longitude
    long1 = long0   dist_x / R / Math.PI * 180

    // Calculate the new latitude
    lat1 = lat0   dist_y / R / Math.PI * 180

    return [long1, lat1]
}