Calculating distance travelled along a flight trajectory
LaMP
Mon, 08/25/2014 - 09:39 am
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.
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
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
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]andVariable LON_ORIGIN = [value], change the LAT[i] and LON[i] to those origin variables, changei+1toi, 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
Yes, it works very well. Thank you very much! I really appreciate that.
August 27, 2014 at 06:53 am - Permalink
I think that in this relation Distance[i] = 2*atan2(sqrt(Distance[i]),sqrt(1-Distance[i])) we need to multiply by the earth radius as well (as it is mentioned here http://en.wikipedia.org/wiki/Haversine_formula).
June 9, 2015 at 02:37 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
// 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
// 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.
// 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