4

I want to calculate the distance from a point given by latitude and longitude to a line-segment given by two points (part of a great-circle). All coordinates are given in WGS84.

enter image description here

I know how to calculate this in Cartesian coordinates but not on a sphere. Can somebody please provide the formula?

Spektre
  • 49,595
  • 11
  • 110
  • 380
user2033412
  • 1,950
  • 2
  • 27
  • 47

2 Answers2

3

  1. spherical to 2D Cartesian

    if the distance is not too far and not around pole singularities you can convert both line segment and line emitting from your point and perpendicular to your segment to spherical coordinates (if they are not already) and use the 2 angles as Cartesian space (ignoring radius).

    1. compute intersection point

    2. convert back to spherical

    3. compute arclength between point and intersection

      hard to say if you're using sphere or WGS84 or what ....

  2. Cartesian 3D

    lets have arc segment AB, sphere of radius R and center C (ideally (0,0,0)) and point P then I see it like this:

    3D cartesian

    1. find intersection point P' between plane ABC and its normal going through point P in 3D Cartesian

    2. project it back on sphere surface

      For spherical surface is this easy as the projection means just to change the vector P'C length to R (if the sphere is centered around (0,0,0)).

      P'' = (R*(P'-C)/|P'-C|) + C

    3. compute arclength between the 2 points |P-P''|

      Also simple for spherical surface just compute the angle between vectors P-C and P''-C

      ang = acos(dot(P-C,P''-C)/(R*R)); // [radians]

      and convert to arclength

      d = ang*R; // [same Units as R]
Spektre
  • 49,595
  • 11
  • 110
  • 380
  • I'm using WGS84. – user2033412 Aug 09 '17 at 10:49
  • @user2033412 that is quite a bit more complicated ... The projection is done iteratively to increase precision ... you simply search `lat` until it fits `x,y,z` you start with spherical projection and then change `lat` to minimize distance ... – Spektre Aug 09 '17 at 10:55
  • @user2033412 see [How to convert a spherical velocity coordinates into cartesian](https://stackoverflow.com/a/41161714/2521214) and the linked answers for some additional ideas ... – Spektre Aug 09 '17 at 11:08
  • 1
    I'm having trouble getting this to work. Take a simple case, see visualization here: https://www.desmos.com/calculator/7vknrruqwc : segment endpoint A is at (1,0,0), segment endpoint B is at (0,1,0), and the point P is at (4/5, 3/5, 0). All are coplanar, so it reduces to a circle, we can see that P is on arc segment AB, so d should be 0. But let's work it out using your method (as I read it). The shortest closest point on straight line AB to P is point I: (3/5,2/5,0). Extend this out radially and you get point P': (2/√13, 3/√13). This is not the same as the original P, so d is nonzero! – hypehuman Sep 02 '20 at 10:01
  • @hypehuman +1 hmm youre right I had bug in the first step ... the intersection must be done between the normal and plane instead of line ... I repaired the answer and updated image. In your case points `P,P',P''` are identical (so the distance is zero) as the `P` already lies on the arc `AB`. the `P'` does no need to lie on the AB line ... that was the bug – Spektre Sep 28 '20 at 09:47
2

This is cross-track distance described here

dxt = asin( sin(δ13) ⋅ sin(θ13−θ12) ) ⋅ R
    where
 δ13 is (angular) distance from start point to third point
     θ13 is (initial) bearing from start point to third point
     θ12 is (initial) bearing from start point to end point
     R is the earth’s radius

You can calculate needed distance and bearings using formulas from given page

distance
a = sin²(Δφ/2) + cos φ1 ⋅ cos φ2 ⋅ sin²(Δλ/2)
c = 2 ⋅ atan2( √a, √(1−a) )
d = R ⋅ c
where   
 φ is latitude, λ is longitude, R is earth’s radius (mean radius = 6,371km);

bearing
θ = atan2( sin Δλ ⋅ cos φ2 , cos φ1 ⋅ sin φ2 − sin φ1 ⋅ cos φ2 ⋅ cos Δλ )
    where
φ1,λ1 is the start point, φ2,λ2 the end point (Δλ is the difference in longitude)

note that angles need to be in radians to pass to trig functions

MBo
  • 77,366
  • 5
  • 53
  • 86