# Projecting Cross-Track distance on great Circle?

I have a route formed by some waypoints (black path) and I would like to project it over the great circle (blue path) connecting the starting and ending waypoint. With the work "project" I mean to find the Latitude and Longitude of the red points.

The red lines are the shortest distances between the black points and the great circle and I am able to calculate it as suggested by this website (section cross-track distance).

Can you tell me how I can calculate the geographical coordinates of the red points, if the only information provided is about the black route and the great circle? • What does your code so far look like?
– PolyGeo
Sep 5, 2016 at 9:21
• cannot provide the code due to IP reasons; I can tell you tho that I have a list of waypoints (black points) and I am able to calculate the red distance by using the formulas provided by the website I linked.., The along-track distance, from the start point to the closest point on the path to the third point, is Formula: dat = acos( cos(δ13) / cos(δxt) ) ⋅ R, where: δ13 is (angular) distance from start point to third point, δxt is (angular) cross-track distance, R is the earth’s radius. Sep 5, 2016 at 9:33

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.

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,result1,lat3,lon3,r=1)
d2=greatCircle(result2,result2,lat3,lon3,r=1)

if d1>d2:
return result2,result2
else:
return result1,result1

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:

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

'''

if r is None:

# 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,p1,p2,p2,p3,p3)
print 'pp',pp

dxt=getCrossTrackDistance(p1,p1,p2,p2,p3,p3,r=1)
print 'dxt',dxt

dat=getAlongTrackDistance(p1,p1,p2,p2,p3,p3,r=1)
print 'dat',dat

dxt2=greatCircle(pp,pp,p3,p3,r=1)
print 'dxt2',dxt2

dat2=greatCircle(pp,pp,p1,p1,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.

• @FedericoGentile I think my answer is wrong. Let me think it through and come back later. Nov 27, 2018 at 9:05
• Ok good to know!! Nov 27, 2018 at 13:10
• @Jason where is the implementation of getAlongTrackDistance ? Can you please share that code? also is your implementation is any different from alongTrackDistanceTo in this page movable-type.co.uk/scripts/latlong.html
– HaBo
Feb 21, 2020 at 15:20
• @HaBo I literally just translated the equation listed under the "Cross-track distance" heading of the https://www.movable-type.co.uk/scripts/latlong.html into python code, so there are no differences. Feb 22, 2020 at 5:24

With PyQGIS it is not difficult to do that. For example, for this situation: where orange line represents a great circle, next code find these "red points" by using 'closestSegmentWithContext' of QgsGeometry class.

``````mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

#for point layer
feats_points = [  feat for feat in layers.getFeatures()  ]

#for great circle
feat_gc = layers.getFeatures().next()

lines  = []
red_points = []

for feat in feats_points:
lines.append([feat.geometry().asPoint(), feat_gc.geometry().closestSegmentWithContext(feat.geometry().asPoint())])
red_points.append(feat_gc.geometry().closestSegmentWithContext(feat.geometry().asPoint()))

print QgsGeometry.fromMultiPoint(red_points).exportToWkt()
print QgsGeometry.fromMultiPolyline(lines).exportToWkt()
``````

After running the above code at the Python Console of QGIS, you can observe printed there the coordinates of "red points" and closest segments (from track green points); both in WKT format.

These layers were displayed with the help of QuickWKT plugin at above image.

I found the answer to my question by using the algorithms suggested on this website.

The idea is to find the along-cross distance defined as:

along-track distance: distance between the start point and the closest point on the path to the third point

In simple words is the length of the great circle between the starting point (blue dot in my drwaing) and the generic projection along it (red dot in my drawing).

Since the goal is to find the coordinates (lat/lon) of the generic projected point along the great circle, we need to combine the following information:

• Lat/lon of the starting point
• Initial bearing angle
• along-track distance

Such inputs allow to calculate the geographical coordinates of the projection along the great circle of a generically located point.

All the equations and algorithms are explained in the website previously linked.