EDIT: I deleted my previous answer as it was wrong. First big mistake I made was performing dot and cross products using spherical coordinates. One needs to convert them to Cartesian first.
(I didn't figure out how to type in math, so I'm pasting images showing all the math).
Below is the python implementation:
def spherical2Cart(lat,lon):
clat=(90-lat)*np.pi/180.
lon=lon*np.pi/180.
x=np.cos(lon)*np.sin(clat)
y=np.sin(lon)*np.sin(clat)
z=np.cos(clat)
return np.array([x,y,z])
def cart2Spherical(x,y,z):
r=np.sqrt(x**2+y**2+z**2)
clat=np.arccos(z/r)/np.pi*180
lat=90.-clat
lon=np.arctan2(y,x)/np.pi*180
lon=(lon+360)%360
return np.array([lat,lon,np.ones(lat.shape)])
def greatCircle(lat1,lon1,lat2,lon2,r=None,verbose=False):
'''Compute the great circle distance on a sphere
<lat1>, <lat2>: scalar float or nd-array, latitudes in degree for
location 1 and 2.
<lon1>, <lon2>: scalar float or nd-array, longitudes in degree for
location 1 and 2.
<r>: scalar float, spherical radius.
Return <arc>: great circle distance on sphere.
'''
if r is None:
r=6371 # km
d2r=lambda x:x*np.pi/180
lat1,lon1,lat2,lon2=map(d2r,[lat1,lon1,lat2,lon2])
dlon=abs(lon1-lon2)
numerator=(cos(lat2)*sin(dlon))**2 + \
(cos(lat1)*sin(lat2) - sin(lat1)*cos(lat2)*cos(dlon))**2
numerator=np.sqrt(numerator)
denominator=sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(dlon)
dsigma=np.arctan2(numerator,denominator)
arc=r*dsigma
return arc
def getCrossTrackPoint(lat1,lon1,lat2,lon2,lat3,lon3):
'''Get the closest point on great circle path to the 3rd point
<lat1>, <lon1>: scalar float or nd-array, latitudes and longitudes in
degree, start point of the great circle.
<lat2>, <lon2>: scalar float or nd-array, latitudes and longitudes in
degree, end point of the great circle.
<lat3>, <lon3>: scalar float or nd-array, latitudes and longitudes in
degree, a point away from the great circle.
Return <latp>, <lonp>: latitude and longitude of point P on the great
circle that connects P1, P2, and is closest
to point P3.
'''
x1,y1,z1=spherical2Cart(lat1,lon1)
x2,y2,z2=spherical2Cart(lat2,lon2)
x3,y3,z3=spherical2Cart(lat3,lon3)
D,E,F=np.cross([x1,y1,z1],[x2,y2,z2])
a=E*z3-F*y3
b=F*x3-D*z3
c=D*y3-E*x3
f=c*E-b*F
g=a*F-c*D
h=b*D-a*E
tt=np.sqrt(f**2+g**2+h**2)
xp=f/tt
yp=g/tt
zp=h/tt
result1=cart2Spherical(xp,yp,zp)
result2=cart2Spherical(-xp,-yp,-zp)
d1=greatCircle(result1[0],result1[1],lat3,lon3,r=1)
d2=greatCircle(result2[0],result2[1],lat3,lon3,r=1)
if d1>d2:
return result2[0],result2[1]
else:
return result1[0],result1[1]
def getCrossTrackDistance(lat1,lon1,lat2,lon2,lat3,lon3,r=None):
'''Compute cross-track distance
<lat1>, <lon1>: scalar float or nd-array, latitudes and longitudes in
degree, start point of the great circle.
<lat2>, <lon2>: scalar float or nd-array, latitudes and longitudes in
degree, end point of the great circle.
<lat3>, <lon3>: scalar float or nd-array, latitudes and longitudes in
degree, a point away from the great circle.
Return <dxt>: great cicle distance between point P3 to the closest point
on great circle that connects P1 and P2.
NOTE that the sign of dxt tells which side of the 3rd point
P3 is on.
'''
if r is None:
r=CONS.EARTH_RADIUS
# get angular distance between P1 and P3
delta13=greatCircle(lat1,lon1,lat3,lon3,r=1.)
# bearing between P1, P3
theta13=getBearing(lat1,lon1,lat3,lon3)*np.pi/180
# bearing between P1, P2
theta12=getBearing(lat1,lon1,lat2,lon2)*np.pi/180
dtheta=np.arcsin(sin(delta13)*sin(theta13-theta12))
dxt=r*dtheta
return dxt
def getAlongTrackDistance(lat1,lon1,lat2,lon2,lat3,lon3,r=None):
'''Compute the distance from the start point to closest point to 3rd point
<lat1>, <lon1>: scalar float or nd-array, latitudes and longitudes in
degree, start point of the great circle.
<lat2>, <lon2>: scalar float or nd-array, latitudes and longitudes in
degree, end point of the great circle.
<lat3>, <lon3>: scalar float or nd-array, latitudes and longitudes in
degree, a point away from the great circle.
Return <dat>: distance from the start point to the closest point on
the great circle connecting P1 and P2.
See also getCrossTrackDistance(), getCrossTrackPoint().
'''
if r is None:
r=CONS.EARTH_RADIUS
# angular distance from P1 to P3
delta13=greatCircle(lat1,lon1,lat3,lon3,r=1.)
# angular distance from Pcloset to P3
dxt=getCrossTrackDistance(lat1,lon1,lat2,lon2,lat3,lon3,r=1)
dat=r*np.arccos(cos(delta13)/cos(dxt/r))
return dat
And here is a test:
p1=[30,100]
p2=[50,210]
p3=[40,180]
pp=getCrossTrackPoint(p1[0],p1[1],p2[0],p2[1],p3[0],p3[1])
print 'pp',pp
dxt=getCrossTrackDistance(p1[0],p1[1],p2[0],p2[1],p3[0],p3[1],r=1)
print 'dxt',dxt
dat=getAlongTrackDistance(p1[0],p1[1],p2[0],p2[1],p3[0],p3[1],r=1)
print 'dat',dat
dxt2=greatCircle(pp[0],pp[1],p3[0],p3[1],r=1)
print 'dxt2',dxt2
dat2=greatCircle(pp[0],pp[1],p1[0],p1[1],r=1)
print 'dat2',dat2
where getCrossTrackDistance() is using the same method as mentioned by OP here, same as getAlongTrackDistance(). The results shown:
pp (56.66932839386002, 185.31664265798963)
dxt 0.297182506587
dat 1.09661554384
dxt2 0.297182506587
dat2 1.09661554384
so we have consistent results.
Hopefully I got it right this time.
EDIT: uploaded getCrossTrackDistance() and getAlongTrackDistance() function definitions.
Best Answer
You can use GeographicLib, which has bindings to several programming languages. For example, the JavaScript version has an example to compute waypoints between two points.
Essentially, you create an inverse geodesic line between the two locations, then generate as many points (either total number or interval distance) between the two points along the geodesic line.