# 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 '16 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. – Federico Gentile Sep 5 '16 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.

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

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

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.

• @FedericoGentile I think my answer is wrong. Let me think it through and come back later. – Jason Nov 27 '18 at 9:05
• Ok good to know!! – Federico Gentile Nov 27 '18 at 13:10

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.