2

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?

enter image description here

  • 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
3

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

enter image description here enter image description here

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]

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.

  • @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
0

With PyQGIS it is not difficult to do that. For example, for this situation:

enter image description here

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.

enter image description here

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

0

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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