Calculating distance travelled along a flight trajectory

Share icon

Hello,
I have a series of longitude, latitude, and altitude points corresponding to a flight trajectory (attached). I would like to calculate the distance (m/km) travelled at each point along the flight trajectory from a fixed point (e.g. the center of Paris 48.85 N, 2.34E). Is this possible to do this with IGOR?
Thank you in advance. Example_flighttrajectory.pxp (9.58 KB)
I'm not at my computer at the moment, but you should check out Vincenty's formulae. Do an Internet search. It shows how to calculate Great Circles between points, which you need in order to properly calculate the distance between two points. I don't think Igor has a built in function to do that, but I could be wrong.

That in mind, I don't understand the part of your question, "from a fixed point." Wouldn't you just want to calculate the distance between each point, and then add them all up to get the total distance? Or are you saying you want the distance from an a priori point to each of the points in your track?

August 25, 2014 at 05:45 pm - Permalink

Hello,
Thank you for your comment and suggestions.

""That in mind, I don't understand the part of your question, "from a fixed point." Wouldn't you just want to calculate the distance between each point, and then add them all up to get the total distance? Or are you saying you want the distance from an a priori point to each of the points in your track?""

I would like to calculate the distance that the plane has travelled from the centre of Paris, which could be calculated by adding the distances between each point and add them together.

Someone had given me a bit of IGOR code a few years ago that did this, but I unfortunately cannot find it. I do not think that I could write this code myself.
I will keep searching...

August 26, 2014 at 06:30 am - Permalink

Does this help? It's code I wrote that I'm pretty sure uses Vincenty's formulae. The LAT is the latitude wave, LON is longitude wave. It calculates the distances between each successive point, storing them to "dist_calc" and "dist_m". It's not elegant, but it does the job.

Macro CalculateDistances()
 
Silent(1)
PauseUpdate
 
Variable i
Make/O/N=(numpnts(LAT)-1) dist_calc; dist_calc = NaN
do
    dist_calc[i] = (sin((LAT[i]-LAT[i+1])*pi/180/2))^2 + cos(LAT[i]*pi/180)*cos(LAT[i+1]*pi/180)*(sin((LON[i]-LON[i+1])*pi/180/2))^2
    dist_calc[i] = 2*atan2(sqrt(dist_calc[i]),sqrt(1-dist_calc[i]))
    i += 1
while(i < numpnts(LAT)-1)
 
//Save as both meters and feet.
Duplicate/O dist_calc dist_m
dist_m = dist_calc * 6371*1000  //convert to meters
dist_calc *= 3959*5280          //convert to feet
 
//Integrate along the path, such that the last point in dist_m_INT will be the total length travelled in meters.
Duplicate/O dist_m dist_m_INT
Integrate dist_m/D=dist_m_INT
 
End


If you wanted to modify this to do it from a fixed point each time, then you would set a Variable LAT_ORIGIN = [value] and Variable LON_ORIGIN = [value], change the LAT[i] and LON[i] to those origin variables, change i+1 to i, and run it. You'd want to get rid of the do-while loop in that case, using Igor's matrix notation, but again, without that it'd still get the job done.

August 26, 2014 at 04:20 pm - Permalink

Hello,
Yes, it works very well. Thank you very much! I really appreciate that.

astrostu wrote:
Does this help? It's code I wrote that I'm pretty sure uses Vincenty's formulae. The LAT is the latitude wave, LON is longitude wave. It calculates the distances between each successive point, storing them to "dist_calc" and "dist_m". It's not elegant, but it does the job.

*****IGOR CODE****
If you wanted to modify this to do it from a fixed point each time, then you would set a Variable LAT_ORIGIN = [value] and Variable LON_ORIGIN = [value], change the LAT[i] and LON[i] to those origin variables, change i+1 to i, and run it. You'd want to get rid of the do-while loop in that case, using Igor's matrix notation, but again, without that it'd still get the job done.

August 27, 2014 at 06:53 am - Permalink

Hi

While searching for a solution to this myself I came across this thread. I'll add my own solutions here for completeness if anyone has need for it. 

Firstly, a function that takes lat,lon, plus optional start and end points, and altitude. The total distance can be calculated in 3 ways:

1. total distance over WGS-84 ellipsoid

2. total distance over WGS-84 ellipsoid accounting for altitude using simple pythagorean theorem - beware this works OK if you have high time resolution positional data, but not necessarily if you have two points far apart.

3. total spherical distance

 

// Calculate total distance along a set of lat-lon co-ordinates
// method=0 (default) - total distance over WGS-84 ellipsoid accounting for altitude above if given
// method=1 - total spherical distance 
Function DistanceAlongTrack(lat,lon[,stPnt,endPnt,method,alt])
    Wave lat,lon
    Variable stPnt,endPnt,method
    Wave alt
    
    // if no altitude is given make a zero wave in its place
    if (paramIsDefault(alt))
        Make/O/FREE/N=(dimSize(lat,0)) alt=0
        endif
    
    // Lat and Lon waves should be the same size 
    if (dimSize(lat,0)!=dimSize(lon,0))
        return 0
        endif
    
    // if no endpoint is given  
    if (paramisDefault(endPnt))
        endPnt=dimSize(lat,0)
        endif
 
    Variable i,totDist
    for (i=stPnt;i<(endPnt-1);i+=1)
        if (method==0) // total distance over WGS-84 ellipsoid accounting for altitude above
            totDist+=EllipsoidalDistanceBetweenTwoPoints(lat[i],lon[i],lat[i+1],lon[i+1],alt1=alt[i],alt2=alt[i+1])
        elseif (method==1) // total spherical distance      
            totDist+=SphericalDistanceBetweenTwoPoints(lat[i],lon[i],lat[i+1],lon[i+1])
        endif
    endfor
    
    // distance in meters   
    return totDist
end

 Secondly there is code to calculate total spherical distance using the Haversine formula

// calculates spherical distance between two points on earth given coords in decimal degrees
// uses haversine formula
// radius of earth is taken to be 6371 km (IUGG/WGS-84 mean)
Function SphericalDistanceBetweenTwoPoints(lat1, long1, lat2, long2)
    Variable lat1, long1, lat2, long2
    
    // convert inputs in degrees to radians
    lat1*= pi/180
    long1*= pi/180
    lat2*= pi/180
    long2*= pi/180
    
    Variable dlon = long1 - long2
    Variable dlat = lat1 - lat2
    Variable aa = ((sin(dlat / 2) ^ 2) + cos(lat1) * cos(lat2) * (sin(dlon / 2) ^ 2))
    Variable cc = 2 * atan2(sqrt(aa), sqrt(1 - aa))
    Variable dd = 6371 * cc
    
    // return distance in metres
    return dd * 1000
end

Thirdly there is code to calculate the distance over a ellipsoid using Vincenty's formula, with optional complication of accounting for altitude.

// calculates ellipsoidal distance between two points on earth given coords in decimal degrees
// uses Vincenty's formula + optional approximately correction for altitude above the ellipsoid 
Function EllipsoidalDistanceBetweenTwoPoints(lat1, long1, lat2, long2 [,alt1, alt2])
    Variable lat1, long1, lat2, long2, alt1, alt2
    
    // convert inputs in degrees to radians
    lat1*= pi/180
    long1*= pi/180
    lat2*= pi/180
    long2*= pi/180
    
    // approximately correct for altitude above the ellipsoid 
    variable delta_alt = abs(alt1 - alt2)
    variable avg_alt = (alt1 + alt2) / 2
    
    // correct for errors at exact poles by adjusting 0.6 millimeters
    if (abs(pi/2-abs(lat1)) < 1e-10)
        lat1 = sign(lat1)*(pi/2-(1e-10))
    endif
    if (abs(pi/2-abs(lat2)) < 1e-10)
        lat2 = sign(lat2)*(pi/2-(1e-10))
    endif
    
    Variable a = 6378137.0 // equatorial radius in meters 
    Variable f = 1/298.257223563 // ellipsoid flattening 
    Variable b = (1 - f)*a 
    Variable tolerance = 1e-11 // to stop iteration
    
    Variable phi1 = lat1 
    Variable phi2 =lat2
    Variable U1 = atan((1-f)*tan(phi1))
    Variable U2 = atan((1-f)*tan(phi2))
    Variable L1 = long1
    Variable L2 = long2
    Variable L = L2 - L1
    
    Variable lambda_old = L + 0
    Variable t
    Variable sin_sigma
    Variable cos_sigma
    Variable sigma
    Variable sin_alpha
    Variable cos_sq_alpha
    Variable cos_2sigma_m
    Variable C
    Variable lambda_new
    
    do  
        t = (cos(U2)*sin(lambda_old))^2
        t += (cos(U1)*sin(U2) - sin(U1)*cos(U2)*cos(lambda_old))^2
        sin_sigma = t^0.5
        cos_sigma = sin(U1)*sin(U2) + cos(U1)*cos(U2)*cos(lambda_old)
        sigma = atan2(sin_sigma, cos_sigma) 
        
        sin_alpha = cos(U1)*cos(U2)*sin(lambda_old) / sin_sigma
        cos_sq_alpha = 1 - sin_alpha^2
        cos_2sigma_m = cos_sigma - 2*sin(U1)*sin(U2)/cos_sq_alpha
        C = f*cos_sq_alpha*(4 + f*(4-3*cos_sq_alpha))/16
        
        t = sigma + C*sin_sigma*(cos_2sigma_m + C*cos_sigma*(-1 + 2*cos_2sigma_m^2))
        lambda_new = L + (1 - C)*f*sin_alpha*t
        
         // correct for convergence failure in the case of essentially antipodal points
        if (lambda_new > pi)
            Print "Points are essentially antipodal. Precision may be reduced slightly."
            lambda_new = pi
            break
        endif
            
        if (abs(lambda_new - lambda_old) <= tolerance)
            break
        else
            lambda_old = lambda_new
        endif
    while (1)       
    
    u2 = cos_sq_alpha*((a^2 - b^2)/b^2)
    variable Av = 1 + (u2/16384)*(4096 + u2*(-768+u2*(320 - 175*u2)))
    variable Bv = (u2/1024)*(256 + u2*(-128 + u2*(74 - 47*u2)))
    t = cos_2sigma_m + 0.25*Bv*(cos_sigma*(-1 + 2*cos_2sigma_m^2))
    t -= (Bv/6)*cos_2sigma_m*(-3 + 4*sin_sigma^2)*(-3 + 4*cos_2sigma_m^2)
    Variable delta_sigma = Bv * sin_sigma * t
    Variable s = (b+(avg_alt-(delta_alt/2)))*Av*(sigma - delta_sigma)
    
    // pythagorean theorem
    s = sqrt(s^2 + delta_alt^2)
    
    // in meters    
    return s
end

Some of this is converted from elsewhere and different languages and I take no credit for it.

September 22, 2020 at 06:21 am - Permalink