Projecting Cross-Track distance on great Circle?

Question

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?

Answer 1

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.

Answer 2

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[0].getFeatures()  ] 

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

lines  = []
red_points = []

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

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.

Answer 3

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, it 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.

Answer 4

There is an intersection(p1, brng1, p2, brng2) function on https://www.movable-type.co.uk/scripts/latlong.html, p1 is start, brng1 is the bearing from start to end, p2 is a given black point, brng2 is brng1 +/- 90 based on the direction, and you should be able to get the red dot. Notice that you don't even need to calculate the cross track distance.