-1

I am using the Haversine formula to calculate the distance from two latitude-longitude pairs.

function getDistanceFromLatLonInKm(lat1,lon1,lat2,lon2) {
    var R = 6371; // Radius of the earth in km
    var dLat = deg2rad(lat2-lat1);
    var dLon = deg2rad(lon2-lon1); 
    var lat1 = deg2rad(lat1);
    var lat2 = deg2rad(lat2);

    var a = 
        Math.sin(dLat/2) * Math.sin(dLat/2) + 
        Math.sin(dLon/2) * Math.sin(dLon/2) * 
        Math.cos(lat1) * Math.cos(lat2);

    var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a)); 
    var d = R * c;
    return d;
}

Given a starting point (lat1, lat2), the distance required to move on a straight line and the angle, I need to determine the endpoint (as in lat2 and lon2).

See my attempt below:

function getFinalLatLon(lat1, lon1, distance, angle) {
    var R = 6371; // Radius of the earth in km
    var c = distance/R;
    // Math.atan2(Math.sqrt(a), Math.sqrt(1-a)) = c/2
    var a = // stuck here

    // looking for this part of the code

    return [lat2, lon2];
}
Timothy
  • 4,198
  • 6
  • 49
  • 59

2 Answers2

4

If you are moving horizontally, you can increment the longitude by distance / (R * cos(lat)). No atan needed.


EDIT: Since you wanted a formula for the general case, consider the following geometric derivation:

  • Front view:

    enter image description here

  • Side view:

    enter image description here

  • Entire setup:

    enter image description here

Notes:

  • r is the unit vector of your starting position, and s is the endpoint.

    enter image description here

  • a, b, c are intermediate vectors to aid calculation.

  • (θ, φ) are the (lat, long) coordinates.
  • γ is the bearing of the direction you are going to travel in.
  • δ is the angle travelled through (distance / radius R = 6400000m).

We need a, b to be perpendicular to r and also a aligned with North. This gives:

enter image description here

c is given by (simple trigonometry):

enter image description here

And thus we get s (through some very tedious algebra):

enter image description here

Now we can calculate the final (lat, long) coordinates of s using:

enter image description here


Code:

function deg2rad(deg) { return deg * (Math.PI / 180.0) }
function rad2deg(rad) { return rad * (180.0 / Math.PI) }

function getFinalLatLong(lat1, long1, distance, angle, radius) {
    // calculate angles
    var delta = distance / radius,
        theta = deg2rad(lat1),
        phi   = deg2rad(long1),
        gamma = deg2rad(angle);

    // calculate sines and cosines
    var c_theta = Math.cos(theta), s_theta = Math.sin(theta);
    var c_phi   = Math.cos(phi)  , s_phi   = Math.sin(phi)  ;
    var c_delta = Math.cos(delta), s_delta = Math.sin(delta);
    var c_gamma = Math.cos(gamma), s_gamma = Math.sin(gamma);

    // calculate end vector
    var x = c_delta * c_theta * c_phi - s_delta * (s_theta * c_phi * c_gamma + s_phi * s_gamma);
    var y = c_delta * c_theta * s_phi - s_delta * (s_theta * s_phi * c_gamma - c_phi * s_gamma);
    var z = s_delta * c_theta * c_gamma + c_delta * s_theta;

    // calculate end lat long
    var theta2 = Math.asin(z), phi2 = Math.atan2(y, x);

    return [rad2deg(theta2), rad2deg(phi2)];
}

Test cases:

  1. Input (lat, long) = (45, 0), angle = 0, distance = radius * deg2rad(90) => (45, 180) (as I said before)

  2. Input (lat, long) = (0, 0), angle = 90, distance = radius * deg2rad(90) => (0, 90) (as expected - start at equator, travel east by 90 longitude)

  3. Input (lat, long) = (54, 29), angle = 36, distance = radius * deg2rad(360) => (54, 29) (as expected - start at any random position and go full-circle in any direction)

  4. Interesting case: input (lat, long) = (30, 0), everything else same. => (0, 90) (we expected (30, 90)? - not starting at equator, travel by 90 degrees to North)

    The reason for this is that 90 degrees to North is not East (if you're not at the equator)! This diagram should show why:

    enter image description here

    As you can see, the path of movement at 90 degrees to the North is not in the direction of East.

meowgoesthedog
  • 14,670
  • 4
  • 27
  • 40
  • Let me prove this one before I accept... Sounds so simple it blows my mind. – Timothy Jul 21 '17 at 10:24
  • @TechyTimo oops hang on, sorry, there was a bit missing – meowgoesthedog Jul 21 '17 at 10:26
  • 1
    @meowgoesthedog Multiply by 180/Pi to get longitude increment in degrees – MBo Jul 21 '17 at 10:44
  • Great explanation. Seems to fail on the pass the first test case though - Taking earth radius = 6371000 metres, distance = radius * deg2rad(90) and getFinalLatLong(45, 0, distance, 90, radius) returns [4.96.., 90] – Timothy Jul 22 '17 at 18:02
  • 1
    @TechyTimo you didnt show the rest of that `lat` value though - 4.961562726608714e-15 - i.e. a *very* small number close to the expected value of zero. This is a floating point error typical to these types of calculations, and is approximately the *machine epsilon* of double precision numbers. – meowgoesthedog Jul 22 '17 at 18:22
  • Maybe you can correct this with adding decimal point precision? toFixed(10) maybe? The second test failed too but the third and forth were close. Maybe provide a fiddle with working tests. – Timothy Jul 22 '17 at 18:51
  • @TechyTimo this is not a failure - encountering small errors like this are perfectly normal in floating point calculations. I imagine you got something like .999999999999992 in the tests afterwards - again these values are well with the margin of error – meowgoesthedog Jul 22 '17 at 18:54
  • The values I am getting are way off your projections. This clearly cannot be a small error. Do a JSfiddle to complete your answer. – Timothy Jul 23 '17 at 06:36
0

I just found a similar question here and I followed the solution to come up with a function that works for my case.

Hope it helps someone else:

function getFinalLatLon(lat1, lon1, distance, angle){ 
    function deg2rad(deg) {
        return deg * (Math.PI/180)
    }
    // dy = R*sin(theta) 
    var dy = distance * Math.sin(deg2rad(angle)) 
    var delta_latitude = dy/110574
    // One degree of latitude on the Earth's surface equals (110574 meters
    delta_latitude = parseFloat(delta_latitude.toFixed(6));

    // final latitude = start_latitude + delta_latitude
    var lat2 = lat1 + delta_latitude

    // dx = R*cos(theta) 
    var dx = distance * Math.cos(deg2rad(angle)) 
    // One degree of longitude equals 111321 meters (at the equator)
    var delta_longitude = dx/(111321*Math.cos(deg2rad(lat1))) 
    delta_longitude = parseFloat(delta_longitude.toFixed(6));

    // final longitude = start_longitude + delta_longitude
    var lon2 = lon1 + delta_longitude

    return [lat2, lon2];
}

The angle is 0 degrees for a horizontal move. You can switch that as you wish. If someone is moving north that would be 90 deg. 135 degrees for north west and so on...

Timothy
  • 4,198
  • 6
  • 49
  • 59
  • Not sure if this works. Say if you move due North by `2 * (90 - lat)`. You should end up at `lat` again, but you get `180 - lat`, which is incorrect. – meowgoesthedog Jul 21 '17 at 11:07
  • I think its fixed - I was missing a degrees to radians conversion that was very much required. – Timothy Jul 22 '17 at 04:09
  • although that is necessary, it doesn't solve the issue I pointed out above. – meowgoesthedog Jul 22 '17 at 10:38
  • I dont understand your question.. You move North by what? Can you use actual figures... – Timothy Jul 22 '17 at 10:40
  • say if you are at latitude = 45. If you move North by distance (earth radius) * 90, you should still arrive at latitude 45 (just on the other side of the earth, i.e. your longitude changes by 180). However if you run your code, the latitude becomes 135, and the longitude doesn't change at all; I think this is an incorrect result. – meowgoesthedog Jul 22 '17 at 10:43
  • To sum up my point, I believe the code works fine *until* you cross the coordinate singularity at one of the poles – meowgoesthedog Jul 22 '17 at 10:45
  • In your case, what is "(earth radius) * 90". And why? Please elaborate this. – Timothy Jul 22 '17 at 11:05
  • multiplying the angle you move through (in radians) by the earth's radius gives the distance travelled on the earth's surface; I should have highlighted that 90 is degrees, equal to Pi / 2 radians. Give me a few minutes and I will upload a diagram of the situation – meowgoesthedog Jul 22 '17 at 11:10
  • Let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/149882/discussion-between-techytimo-and-meowgoesthedog). – Timothy Jul 22 '17 at 11:13
  • Hi, I updated my answer with a more complete derivation – meowgoesthedog Jul 22 '17 at 17:41