How to create a point along a line given a distance from another point

Question

I'm trying to figure out a point(D) on B-C that is a specific distance from A.
I really don't have any idea on how to proceed with the ellipsoidal math so I figured I'd ask here. How would I figure out the longitude, latitude of D which is 5mi(26400ft) away from A, on the line B-C? The bearing of A-D would be just as useful as i could figure out the point 5mi away from A with the bearing.

I wouldn't mind if the formula lined up with the Vincenty model of the earth.

Using the Vincenty formula's I've already calculated the following.

Points (longitude, latitude)
A: -86.32899412316736,41.12507719625815
B: -86.29237110757369,41.16845596588609
C: -86.19689279712584,41.16923568742915

A-B Distance: 18750ft
B-C Distance: 26289ft
A-B Initial Bearing: 32.5331889647
B-C Initial Bearing: 89.3493884487

enter image description here

Thanks for the help!

score 1 · Accepted Answer · answered May 25 '13 at 11:54

This is most easily solve using the equidistant azimuthal projection. Guess some point for D, e.g., A. Transform A, B, C, to equidistant azimuthal projection using D as a center. In projected space, solve for D' (i.e., D' lies on BC and is a distance s from A). Project D' back to lat,lon and update D with this position. Repeat. This converges quadratically.

To see this in action, grab http://www.mathworks.com/matlabcentral/fileexchange/39108 and http://www.mathworks.com/matlabcentral/fileexchange/39366 from MathWorks File Exchange.

Write a .m file to solve the planar problem, planesolve.m. Here is a not-very-elegant version

function D = planesolve(pts, s)
% pts is 3x2 matrix of [A;B;C] lying in a plane.  Find point D which lies
% on B-C and which is a distance s from A.  There may be 0, 1, or 2
% solutions.  Here assume there are two and pick the one most in the
% direction B-C.

  xa=pts(1,1); ya=pts(1,2);
  xb=pts(2,1); yb=pts(2,2);
  xc=pts(3,1); yc=pts(3,2);
% write D = t*B + (1-t)*C, solve for t
  t = ((yc-ya)*(yc-yb)+(xc-xa)*(xc-xb) + ...
       sqrt((-xb^2+2*xa*xb-xa^2+s^2)*yc^2+...
            (((2*xb-2*xa)*xc-2*xa*xb+2*xa^2-2*s^2)*yb+...
             ((2*xa-2*xb)*xc+2*xb^2-2*xa*xb)*ya)*yc+...
            (-xc^2+2*xa*xc-xa^2+s^2)*yb^2+...
            (2*xc^2+(-2*xb-2*xa)*xc+2*xa*xb)*ya*yb+...
            (-xc^2+2*xb*xc-xb^2)*ya^2+s^2*xc^2-2*s^2*xb*xc+s^2*xb^2))/...
      ((yc-yb)^2+(xc-xb)^2);
  xd=t*xb+(1-t)*xc; yd=t*yb+(1-t)*yc;
  D=[xd,yd];
end

Then run the following script (in either matlab or octave)

g=[41.12507719625815,-86.32899412316736;...
   41.16845596588609,-86.29237110757369;...
   41.16923568742915,-86.19689279712584];
s = 5*1760*36*0.0254;
gd=g(1,:);
for i=1:5,
  [x,y] = eqdazim_fwd(gd(end,1),gd(end,2),g(:,1),g(:,2));
  d = planesolve([x,y],s);
  [lat,lon] = eqdazim_inv(gd(end,1),gd(end,2),d(1),d(2));
  gd = [gd;lat,lon];
end
format long;
gd

This prints out

gd =

  41.125077196258147 -86.328994123167362
  41.167417462129485 -86.406778034888063
  41.167417401850848 -86.406778092275701
  41.167417401850948 -86.406778092275616
  41.167417401850948 -86.406778092275601
  41.167417401850948 -86.406778092275601

which gives the converging sequence of approximations (lat,lon) for D. (Here the WGS84 ellipsoid is assumed. eqdazi_fwd and eqdazi_inv take an optional argument which lets you specify the ellipsoid.)

This ends up being really close, the distance is correct but It doesn't land on line segment B-C? — Phil, May 26 '13 at 17:41
Yes! The distance is right because the projection is equidistant (with D as the center). D lies on B-C because the projection is azimuthal (with D as the center). I should add that you can pretty much choose any point as a starting point: A, B, C, (B+C)/2, the North Pole, the Eiffel Tower all work fine. — cffk, May 28 '13 at 12:42
the distance is right.. but how would you go about getting a point between B & C. instead of just on the line B-C. — Phil, May 31 '13 at 04:05
My apologies, I got you the solution outside the range B-C. Change the sign of the sqrt in planesolve.m and rerun and you will get the solution you want 41.1687907834032 -86.2525504601801. In general you should modify the code to try both solutions in planesolve and then select the one (if any) which gives you opposite azimuths D->B and D->C. — cffk, May 31 '13 at 09:36

score 0 · Answer 2 · answered May 23 '13 at 20:00

Intersection of two paths given start points and bearings could be used. You would need to check the resulting point to make sure it is between B and C. Since there's no spheroid flattening, I'm not sure how close it would match up with Vincenty's formula.

Formula:    

d12 = 2.asin( √(sin²(Δφ/2) + cos(φ1).cos(φ2).sin²(Δλ/2)) )
φ1 = acos( sin(φ2) − sin(φ1).cos(d12) / sin(d12).cos(φ1) )
φ2 = acos( sin(φ1) − sin(φ2).cos(d12) / sin(d12).cos(φ2) )

if sin(λ2−λ1) > 0
    θ12 = φ1, θ21 = 2.π − φ2
else
    θ12 = 2.π − φ1, θ21 = φ2

α1 = (θ1 − θ12 + π) % 2.π − π
α2 = (θ21 − θ2 + π) % 2.π − π

α3 = acos( −cos(α1).cos(α2) + sin(α1).sin(α2).cos(d12) )
d13 = atan2( sin(d12).sin(α1).sin(α2), cos(α2)+cos(α1).cos(α3) )
φ3 = asin( sin(φ1).cos(d13) + cos(φ1).sin(d13).cos(θ1) )
Δλ13 = atan2( sin(θ1).sin(d13).cos(φ1), cos(d13)−sin(φ1).sin(φ3) )
λ3 = (λ1+Δλ13+π) % 2.π − π
where   

φ1, λ1, θ1 : 1st point & bearing
φ2, λ2, θ2 : 2nd point & bearing
φ3, λ3 : intersection point

% = mod
note –  if sin(α1)=0 and sin(α2)=0: infinite solutions
if sin(α1).sin(α2) < 0: ambiguous solution
this formulation is not always well-conditioned for meridional or equatorial lines

wouldn't that just give point B – Phil May 23 '13 at 20:07 — Phil, May 23 '13 at 20:07

How to create a point along a line given a distance from another point

2 Answers2