Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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!

share|improve this question
add comment

2 Answers

up vote 1 down vote accepted

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.)

share|improve this answer
    
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 4: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 9:36
    
that worked perfectly... thanks! –  Phil Jun 3 '13 at 20:09
add comment

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
share|improve this answer
1  
wouldn't that just give point B –  Phil May 23 '13 at 20:07
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.