# frosch03.de/posts/2015-10-11-little_vor_navigation_helper

## Idea

So this time i had the problem of getting the correct vor radial for navigation. To get from one airport to another, i usually set up a route i want to fly by. My thought was, to create a little program that tells you for each point of the route the two nearest vor stations with the radial. With this information it is then possible to triangulate the point of the route by radio navigation.

A quick search revealed openaip.net, an free list of navigation aids that you could download. You get a download for navigation aids per country. The date is presented as .aip xml file.

I then used .gpx as input format for the flight route, since this format is also xml based.

As it is pretty straigt forward within python to work with data in xml format i choose python as the language for this task. I also use the numpy library as well a few more librarys, so my import section looks like this:

{% highlight python %} import sys import numpy as np import xml.etree.ElementTree as ET from itertools import chain {% endhighlight %}

The task itself is basically not more then calculating for each point of the route the distance to every vor station, sort the result by distance and drop all but the two nearest ones. Calculate the bearing of thouse two also.

A finished version with a little more functionality to use it from a command line is available at github

The navigation aid files could be downloaded from openaip.net. Note that you need to create an account on that page to actually get to the download link.

## Geometry

### Distance

To calculate the distince of two points on an euclidian surface, one has just to follow the theorem of pythagoras. Calculate square root of the sum of the square of the absolute value of the x coordinate differences and the y coordinate differences. Or for short: \[d = \sqrt{|x_1 - x_2|^2 + |y_1 - y_2|^2}\]

Unfortunatly this does not apply on a sphere like the earth. Since i didn't want to derive that formula on my own, i just look it up. To calculeate the distance between two points on a sphere the haversine formula could be used: \[a = \sin\^{2}(\frac{\Delta\phi}{2}) + \cos(\phi/{1}) * \cos(\phi/{2}) * \sin\^{2}(\frac{\Delta\lambda}{2}) \]

\[c = 2 * \arctan2(\sqrt(a), \sqrt(1-a)) \]

\[d = c * R \]

- \(R\) is just the mean radius of the earth: ~6.371km
- \(\phi\) is the value of latitude
- \(\lambda\) is the value of longitude

I belive that the intermediate value \(c\) is used, since it is unit free. Only the last step, by multiplying it with the earth's radius introduces the unit. If you multiply it with the radius measured in metres you end up with the distance in meter. The same is true if you multiply \(c\) with the radius in miles.

All right, so the corresponding python code, that calculates the distance looks pretty similar:

{% highlight python lineno %} def distance((lat1, lon1), (lat2, lon2)):
r\_m = 6.371e6 r\_nm = 3443.92 phi1 = np.radians(lat1) phi2 =
np.radians(lat2) dphi = np.radians(lat2-lat1) dlbd =
np.radians(lon2-lon1) a = np.sin(dphi/2) * np.sin(dphi/2) +

np.cos(phi1) * np.cos(phi2) *

np.sin(dlbd/2) * np.sin(dlbd/2) c = 2 * np.arctan2(np.sqrt(a),
np.sqrt(1-a)) d = r\_nm * c return d {% endhighlight %}

Lines 2 and 3 are just the radius in meters and nautical miles. Within Lines 4-7 i calculate the radians for (φ\_{1,2}\), \(\Delta\phi\) and \(\Delta\lambda\). Line 8-10 calculates then \(a\), \(c\) and \(d\) according to the haversine formulas.

### Bearing

For the bearing, there is the same situation. The euclidian formulas do not apply, so i needed to look up the formula for a sphere.

\[ x\_1 = \sin(\Delta\lambda) * \cos(\phi2)\] \[ x\_2 = \cos(\phi\_1) * \sin(\phi\_2) - \sin(\phi\_1) * \cos(\phi\_2) * \cos(\Delta\lambda)\] \[ brg = \arctan2(x\_1, x\_2) \]

Within the code i also need to first calculate the radian values. Then i calculate the two opperands of the \(\arctan2()\) function and finaly get the result back in degrees. To avoid negative bearings the result is added by 360 and is then divided by flotingpoint modulo 360.

{% highlight python lineno %} def bearing((lat1, lon1), (lat2, lon2)):
phi1 = np.radians(lat1) phi2 = np.radians(lat2) dlbd =
np.radians(lon2-lon1) x1 = np.sin(dlbd)/np.cos(phi2) x2 =
(np.cos(phi1)/np.sin(phi2)) -

(np.sin(phi1)/np.cos(phi2)/np.cos(dlbd)) brg = np.degrees(np.arctan2(x1,
x2)) return np.fmod((360 + brg), 360) {% endhighlight %}

So with these two functions, it is possible to calculate the distance of two points and their relative bearing to each other. From a mathematical standpoint the problem is solved at this point ;-)

## XML Handling

So the basic xml parsing is done by the `xml.etree.ElementTree`

library.
One parses a xml file with `ET.parse('filename.xml')`

and gets back the
root tree from which the root is returned by `getroot()`

. So the navaids
of an .aip file are available via:

{% highlight python lineno %} element\_tree = ET.parse('filename') element\_root = element\_tree.getroot() {% endhighlight %}

### aip files

An .aip file consists of a list of navaids. A navaid has a type and various other elements. It looks like the following:

{% highlight xml lineno %} DE BAYREUTH BAY 49.9866676331 11.6383333206 487 110.06 2.47592 FALSE {% endhighlight %}

Now there are different types of vor's together with ndb typed navaids
within a .aip file. Since i would like to navigate with vors, i need to
fetch all vor's from the .aip file (e.g.
`openaip_navaid_germany_de.aip`

).

{% highlight python lineno %} navaids\_tree = ET.parse('openaip\_navaid\_germany\_de.aip') navaids\_root = navaids\_tree.getroot() {% endhighlight %}

Within an aip file there are vor's of the type: - VOR - DVOR - DVOR-DME

- VORTAC

Therefore fetching all vor's is fetching all four vor types and unioning them together:

{% highlight python lineno %} vor\_navaids =
navaids\_root[0].findall(".//[@TYPE='VOR']")

- navaids\_root[0].findall(".//[@TYPE='DVOR']")
- navaids\_root[0].findall(".//[@TYPE='DVOR-DME']")
- navaids\_root[0].findall(".//[@TYPE='VORTAC']") {% endhighlight %}

The main data type for navigation aids used later on consists of bearing, distance, frequency and id. So i wrote getter functions, to read these information out of aip files.

{% highlight python lineno %} def aipGetLatLon(navaid): lat = float(navaid.find('GEOLOCATION/LAT').text) lon = float(navaid.find('GEOLOCATION/LON').text) return (lat, lon)

def aipGetId(navaid): return (navaid.find('ID').text)

def aipGetFrequency(navaid): return (float(navaid.find('RADIO/FREQUENCY').text)) {% endhighlight %}

With these getters and the list of navaids from the .aip file, the list
of vor's could be calculated in a format, that is much more suitable for
the given task. A tuple of a *description* and the latitude/longitude
will do the job perfectly. The *description* itself is a tuple of an id
string and a frequency float.

((id_string, frequency) , (latitude, longitude))

The following maps a lambda expression over the navaid list and greps the data in exactly this format:

{% highlight python lineno %} vors = map(lambda n: ((aipGetId(n), aipGetFrequency(n)), aipGetLatLon(n)), vor\_navaids) {% endhighlight %}

### gpx files

The flight path or route is stored within an .gpx file. The xml is pretty basic and the main information is stored within a list of route points:

{% highlight xml lineno %} 400.2 A\_EDSP POLTRINGEN {% endhighlight %}

The file is processed just as the .aip file. First the tree parsed and the root returned:

{% highlight python lineno %} route\_tree = ET.parse('route.gpx') route\_root = route\_tree.getroot() {% endhighlight %}

To get the list of route points only the *rtept*'s of the gpx must be
grepped:

{% highlight python lineno %} route = route\_root[0].findall("{http://www.topografix.com/GPX/1/1}rtept") {% endhighlight %}

Also getters for the information of the gpx file are needed:

{% highlight python lineno %} def gpxGetLatLon(pt): lat = float(pt.attrib['lat']) lon = float(pt.attrib['lon']) return (lat, lon)

def gpxGetName(pt): return (pt.findall("{http://www.topografix.com/GPX/1/1}name")[0].text) {% endhighlight %}

The route is then stored internaly in a tuple of an id string and the latitude/longitude information:

(id_string , (latitude, longitude))

This list of route points is created by mapping the corresponding lambda expression over the list of route points:

{% highlight python lineno %} route = map(lambda pt: (gpxGetName(pt), gpxGetLatLon(pt)), route\_) {% endhighlight %}

## Calculations along the Route

### Navigation Aids along the Route

The previous sections showed, how the information of navaids and route
points are accessed. This part shows, how to combine them together to
get for every point on the road, the two nearest vor's with their
distance and their bearing. Since we do not know the nearest vor of each
waypoint, we consider every known vor as *along the route*.

{% highlight python lineno %} def getNavaidsEnroute(navaids, route): return map(lambda (sr, r): map(lambda ((sn, fn), n): ((sr + "->" + sn, fn), calcBrgDst(r, n)), navaids), route) {% endhighlight %}

The function `calcBrgDst`

calculates just the bearing and the distance
of two points and returns both as tuple:

{% highlight python lineno %} def calcBrgDst((lat1, lon1), (lat2, lon2)): brng = bearing((lat1, lon1), (lat2, lon2)) dist = distance((lat1, lon1), (lat2, lon2)) return (brng, dist) {% endhighlight %}

The `getNavaidsEnroute`

function takes a list of navaids and a route. It
then maps, for each route point, the `calcBrgDst`

function with the
route point over each navaid position.

The result is a list of waypoints, each contianing a list of navaid discription together with bearing and distance.

### Uniqify, Sort and Shorten

#### Unique

A sucessful triangulation could only be done with two distinct navigation aids. Therefore it is importend that no navaid is multiple times within the list of navaids along the route. A fast way to unique a list of elements, i found at stackoverflow:

{% highlight python lineno %} def uniqFast(seq): seen = set() seen\_add = seen.add return [x for x in seq if not (x in seen or seen\_add(x))] {% endhighlight %}

#### Sort

A list of everything could be sorted, if one could present a key
function to the library function `sorted`

. All possible key functions
for the navaid data are listed here:

{% highlight python lineno %} def getDistance(brgdst): return brgdst[1][1]

def getBearing(brgdst): return brgdst[1][0]

def getFrequency(brgdst): return brgdst[0][1] {% endhighlight %}

The sort function is then pretty straigt forward:

{% highlight python lineno %} def sortDistance(brgdsts): return sorted(brgdsts, key=getDistance) {% endhighlight %}

#### Shorten

All could be combined to the function `getNearestNavaidsEnroute`

. This
function takes a list of navaids and a list of waypoints and optionally
a count of how many navaids shold be left in the nearest list.

{% highlight python lineno %} def getNearestNavaidsEnroute(navaids, route, count=2): navaids\_enroute = getNavaidsEnroute(navaids, route) return map(lambda naids: (uniqFast(sortDistance(naids)))[:count], navaids\_enroute) {% endhighlight %}

It maps the combined `uniqFast(sortDistance())`

function over the
navaids along the route and drops all but `count`

navaids. This results
in the desired list of waypoints, together with a list of a few nearest
navigation aids.

## Showing the Result

Last but not least, a show function helps to visualize the calculation result:

{% highlight python lineno %} def showNavaidsEnroute(naids, route): for x in naids: for y in x: s = getName(y) b = str(int(getBearing(y))) d = str(round(getDistance(y), 2)) f = "{:7.3f}".format(round(getFrequency(y), 3)) print s + "(" + f + ")"":" + b + u"°," + d + "NM" print "–" {% endhighlight %}

The result of a short trip around stuttgard would result in this output:

A_EDSP->STG(116.850): 53°, 15.32NM A_EDSP->SUL(116.100): 230°, 15.54NM -- B_N 48->SUL(116.100): 254°, 12.78NM B_N 48->STG(116.850): 37°, 19.66NM -- C_N 48->SUL(116.100): 228°, 4.42NM C_N 48->STG(116.850): 52°, 26.43NM -- D_N 48->KRH(115.950): 331°, 10.16NM D_N 48->STG(116.850): 112°, 23.48NM -- E_EDTQ->LBU(109.200): 56°, 5.45NM E_EDTQ->STG(116.850): 172°, 10.1NM --