Oct 11, 2015 - Little VOR Navigation Helper

Comments

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:

import sys
import numpy as np
import xml.etree.ElementTree as ET
from itertools import chain

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:

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

Lines 2 and 3 are just the radius in meters and nautical miles. Within Lines 4-7 i calculate the radians for (\phi_{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.

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)

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:

element_tree = ET.parse('filename')
element_root = element_tree.getroot()

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:

<NAVAID TYPE="VOR">
  <COUNTRY>DE</COUNTRY>
  <NAME>BAYREUTH</NAME>
  <ID>BAY</ID>
  <GEOLOCATION>
    <LAT>49.9866676331</LAT>
    <LON>11.6383333206</LON>
    <ELEV UNIT="M">487</ELEV>
  </GEOLOCATION>
  <RADIO>
    <FREQUENCY>110.06</FREQUENCY>
  </RADIO>
  <PARAMS>
    <DECLINATION>2.47592</DECLINATION>
    <ALIGNEDTOTRUENORTH>FALSE</ALIGNEDTOTRUENORTH>
  </PARAMS>
</NAVAID>

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

navaids_tree = ET.parse('openaip_navaid_germany_de.aip')
navaids_root = navaids_tree.getroot()

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:

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']")

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.

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

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:

vors = map(lambda n: ((aipGetId(n), aipGetFrequency(n)), aipGetLatLon(n)), vor_navaids)

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:

<rtept lat="48.546666666695" lon="8.9450000000061">
    <ele>400.2</ele>
    <name>A_EDSP</name>
    <cmt>POLTRINGEN</cmt>
</rtept>

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

route_tree = ET.parse('route.gpx')
route_root = route_tree.getroot()

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

route_ = route_root[0].findall("{http://www.topografix.com/GPX/1/1}rtept")

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

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)

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:

route = map(lambda pt: (gpxGetName(pt), gpxGetLatLon(pt)), route_)

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

def getNavaidsEnroute(navaids, route):
    return map(lambda (sr, r):
               map(lambda ((sn, fn), n):
                   ((sr + "->" + sn, fn), calcBrgDst(r, n)),
                   navaids),
               route)

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

def calcBrgDst((lat1, lon1), (lat2, lon2)):
    brng = bearing((lat1, lon1), (lat2, lon2))
    dist = distance((lat1, lon1), (lat2, lon2))
    return (brng, dist)

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:

def uniqFast(seq):
    seen = set()
    seen_add = seen.add
    return [x for x in seq if not (x in seen or seen_add(x))]

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:

def getDistance(brgdst):
    return brgdst[1][1]


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


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

The sort function is then pretty straigt forward:

def sortDistance(brgdsts):
    return sorted(brgdsts, key=getDistance)

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.

def getNearestNavaidsEnroute(navaids, route, count=2):
    navaids_enroute = getNavaidsEnroute(navaids, route)
    return map(lambda naids:
               (uniqFast(sortDistance(naids)))[:count],
               navaids_enroute)

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:

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 "--"

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