SimpleFly - GPS Geoid
A flight app needs to be able to provide altitude in Mean Sea Level (MSL). When developing SimpleFly, I found out that the earth is not round! (OK, I kind of knew this already, anyway...) This page describes why the earth is not round and why the sea is not flat. And, if you are interested in MSL, WGS84, and EGM then please keep reading.
But first, if you are interested in the SimpleFly app, the details that can be found here:
https://worktablecnc.us/projects/simplefly.html
GPS Shenanigans
The raw GPS sensor information is in:
- Latitude degrees (0 equator, +90 north pole, -90 south pole).
- Longitude degrees (0 Prime Meridian / England, +180 east, -180 west).
- Altitude in meters (WGS84 ellipsoid datum).
So why is my GPS altitude not MSL?
If you ask someone to draw a picture of the Earth, they draw a circle (a sphere in three-dimensions). But because the Earth is spinning, the centrifugal force makes the equator bulge out. So, A more accurate drawing would be an ellipse (an ellipsoid in three-dimensions). Although this was known in the 19th century, when 20th century rocket scientists got involved they eventually agreed on a specific ellipsoid in 1984. This is the World Geodetic System (WGS) 1984 model, or known as WGS84. The WGS84 is the Earth's ideal ellipsoid. This is what our Global Positioning System (GPS) currently uses as a baseline for Mean Sea Level (MSL) (and mean as in average, not as in not nice).
However, the WGS84 model assumes that the density and thus the gravity of the Earth is homogeneous. It is not. Some parts of the Earth have stronger gravity at the surface than others. And because of that, more sea water gets pulled to those areas. The Mean Sea Level around the Earth is actually is not level, it has gradual hills and valleys, and those are plus or minus 100 meters! So in addition to the WGS84 ideal shape, there are Earth Gravitational Models (EGM) which map these gravitation (and thus MSL) differences. The first was the EGM84 developed with the WGS84, but as surface mapping satellites got better, better models were developed. There is the EGM96, EGM2008, and the EGM2020, each more accurate but more complicated. The odd shape that the EGMs define is termed Geoid.
GPS App Programming
To get a basic moving map with altitude shown in MSL on an Android device, you need the ideal WGS84 altitude provided by the device, and then add the real altitude correction from a EGM Geoid model. The details of this can be found below.
On an iOS device this correction is built in with some I-do-not-know EGM model (as always with iOS devices: hopefully easier but with less control).
WGS and EGM Models
WGS and EGM models and support can be found at the United States Government Office of Geomatics.
A couple of notes about the EGM models and support:
-
The Geoid can be defined by a math equation with a series of numerical coefficients using the spherical harmonics (differential equation calculus), or defined by a simple data set.
-
The EGM84 data set is a text file rounded to the half degree (30 minutes). Each row is a space-separated list of longitude, latitude, and correction. Easy. At times they use scientific notation and the string value contains the letter
E, so a function to convert this file to a list/array needs to watch for that. -
The EGM96 data set is a binary file rounded to the quarter degree (15 minutes). This one seems to be more popular in the examples that I have seen. The correction information is arranged by position in the file; the associated Fortran file contains a header explaining it.
-
Both the EGM96 and EGM84 supplied programs are in Fortran. Readable, but somewhat worthless for anything other than reading the comments.
-
In theory, when using a data set, one should linearly interpolate between two points (if the location is halfway between two points, the correction should be halfway between the two values). In practice, just rounding to the nearest correction will provide about a one meter accuracy. Not accurate enough for surveying, good enough for bouncing around in the sky.
Using the Geoid Files
Although one could use the Geoid files as provided, they are in formats for best practices at the time. It is just as easy to parse the information and reformat it to a JSON file, or whatever you need, so your code is easier to write/read/debug.
Reading the EGM96 Geoid File
The EGM96 data set is binary, and here is a Python file to convert it to text for review:
#!/usr/bin/env python3
with open('WW15MGH.GRD', 'r') as file:
data = str(file.read())
print(data)
On review, the file contains returns and newlines for printing, which is kind of nice for looking at but just one more thing to deal with. It is also a bit confusing to figure out (not impossible, just not efficient for this project). They have grouped the data for printing into 9 sets of 160 values (1440 values == 360 * 4; so a full latitude set grouped strangely), with a 10th set of a single value (WTF?). Weird. By comparison, the EGM84 data set is a no-brainer... so I moved on.
Reading the EGM84 Geoid File
The EGM84 data set is text, but contains the letter E for scientific notation. For this application, converting it to a JSON file would be great! Here is a Python file to convert it:
#!/usr/bin/env python3
# Read the file, convert values from string to floating point.
# Some strings have exponents, such as `0.0000000000E+00` and
# Those need to be fixed:
with open('WWGRID.TXT', 'r') as file:
fileList = []
for line in file:
lineList = []
lineSplit = line.split()
for item in lineSplit:
if len(lineSplit) != 3:
print(f"Error: {line}")
else:
itemSplit = item.split('E')
if len(itemSplit) > 1:
item = float(itemSplit[0]) * ( 10 ** (float(itemSplit[1])) )
else:
item = float(item)
lineList.append(item)
fileList.append(lineList)
# Create a usable file from the list:
# Key = (string) Latitude-Longitude
# Value = (double) Correction
# Longitude is from 0 to 360, change to -180 to 180
# and drop 360; it creates a duplicate
print('{')
for line in fileList:
if line[1] != 360:
longitude = (line[1] - 360) if (line[1] > 180) else line[1]
print("\"{:07.2f} {:07.2f}\": {:},".format(line[0], longitude, line[2]))
print('}')
That will print the formatted output. Then just redirect the output to a file named WWGRID.JSON and you are golden! (OK, only if you delete the comma on the very last entity.)