Freesteel Blog » GPS made bad by greediness

GPS made bad by greediness

Thursday, April 6th, 2017 at 11:22 am Written by:

On my last sorry flight I was carrying a total of four GPS units. One of them is a pile of poo. It’s the bluefly one I programmed to read at the frequency of 10 times a second, so greedy was I for precision and data.


Why do they have that option if it’s a complete waste of my time? All that work learning how to remotely step up its baudrate so that it could transmit data at this higher rate.

Here’s how the error comparisons line up, including the minus 3second correction on the 6030 vario’s timestamps. The “u” value is the number of milliseconds from midnight according to the GPS satellite timestamp.

gnames = [ "gps6030", "gpsxcsoar", "gpsTopShelf", "gpsBlueFly" ]
print("samples, gpsA, gpsB, stddev")
for i in range(len(g)):
    for j in range(i+1, len(g)):
        qA, qB = g[i], g[j]
        dt = 3000 if i == 0 else 0 # observed timeslip on 6030 unit
        dx = qA.x - numpy.interp(qA.u-dt, qB.u, qB.x)
        dy = qA.y - numpy.interp(qA.u-dt, qB.u, qB.y)
        k = (dx**2 + dy**2)[tf0:tf1]
        print(len(k), gnames[i], gnames[j], math.sqrt(k.mean()))
samples gpsA gpsB stddev
400 gps6030 gpsxcsoar 6.086
400 gps6030 gpsTopShelf 7.914
400 gps6030 gpsBlueFly 17.040
550 gpsxcsoar gpsTopShelf 7.845
550 gpsxcsoar gpsBlueFly 19.095
4120 gpsTopShelf gpsBlueFly 17.451

Now I have to program the BlueFly to sample at 200ms like the spare top-shelf GPS does, and see if that fixes it.

Meantime, has is that the accelerometer doing with this not-crap GPS?

We can compute the acceleration from the GPS velocity, which is giving as a 2D direction and velocity.

pV[“t”] = pV.index.asi8/1e9 # needs to be a member column to get the shift function

# GPS velocity record to cartesian
pV[“vx”] = numpy.sin(numpy.radians(pV.deg))*pV.vel
pV[“vy”] = numpy.cos(numpy.radians(pV.deg))*pV.vel
dt = 2

# differenciate to get acceleration
pV[“ax”] = (pV.vx.shift(-dt) – pV.vx.shift(dt)) / ((pV.t.shift(-dt) – pV.t.shift(dt)))
pV[“ay”] = (pV.vy.shift(-dt) – pV.vy.shift(dt)) / ((pV.t.shift(-dt) – pV.t.shift(dt)))

Then we can lightly filter the accelerometer vector (transformed to absolute coordinates by the orientation quaternion) and plot it against this GPS velocity derivation

# filter high frequency noise from acceleration value
fvy = utils.FiltFiltButter(fd.pZ.vy, 0.01)



Given what I’ve had to suffer up till now, I think this is a pretty good match on the data.

As a quick view of the lines of acceleration drawn onto the GPS

fig, ax = plt.subplots(figsize=(15,10))
q = fd.pQ0[t0:t1]  # the GPS top shelf dataframe
plt.plot(q.x, q.y)

# the accelerations pre-transformed by orientation to make them absolute
vx = utils.FiltFiltButter(fd.pZ.vx, 0.01)  
vy = utils.FiltFiltButter(fd.pZ.vy, 0.01)
a00 = InterpT(q, vx)   # interpolate acceleration onto GPS
a01 = InterpT(q, vy)
cc = collections.LineCollection([ ((x,y),(x+l*m0s,y+l*m1s))  for x, y, m0s, m1s in zip(q.x, q.y, a00, a01)], colors="r")


The reason for these little snippets of code is because I know I’m going to lose track of them and forget how to do it, so it will save me time if at the last resort I can fetch it from here.

It wasn’t a long flight, and it didn’t end well due to an encounter with a hedge.


There was a somewhat heavy impact with the ground on the other side of the hedge which broke the glider and seriously dented my pride. I have got to get better at judging this shit. It explains why I am spending the day in Sheffield waiting for someone to correct the damage to one of the 2 most important wires.

Still, the data was amazing:


Consider the barometer as the altitude, where it gives away the fact that I tried to hover the glider (while losing airspeed) before “unexpectedly” losing a huge amount of height on a left turn just before the hedge. No external factors at play other than dimwittedness.

Moving swiftly on, I’ve got hold of the terrain files at around 30m (1arcsecond) squares for this part of the globe and can even plot them:

tile = "N52W004.hgt"  # Or some magic to figure it out from position
ft = open(tile, "rb").read()

x = y = numpy.arange(0, 1+1/3600, 1/3600)
terrX, Y = numpy.meshgrid(x, y)
for i in range(3601):
    for j in range(3601):
        p = (((i) * 3601 + (j)) * 2)  # go to the right spot,
        terrX[i][j] = struct.unpack('>h', ft[p:p+2])[0]  # ">h" is a signed two byte integer


Obvious thing to do is to plot the terrain height under the GPS track with altitude:

pQ = fd.pIGC1
pQ["talt"] = pQ.apply(lambda s: terrX[int((*3600)][int((s.lng - (-4))*3600)], axis=1)
plt.plot(pQ[["alt", "talt"]])


There’s an odd discrepancy at the end where I appear to be underground by 75m, but maybe that’s just bad vertical resolution when you are in a valley. It does show close proximity to the hill shortly before, doing what is known as “scratching” where you fly close to the slope where you think there is remaining lift.

I am keeping my eyes out for any developments I can use of openfoam on terrain to get some insights.

There’s a lot of development of that kind of eddy simulations on terrain for the purpose of siting wind turbines. In the future there could be a realtime simulation of the wind over the terrain flows done somewhere on “the cloud” which a glider could log on to in order to predict where the lift is to get to the clouds.

Some of this sort of lift is large scale, like wave lift, and ought to be manifestly predictable. Likewise, knowing where the zones of sink lie that are known to suppress thermal development.

This is not a big step. We already have the rasp forecasts by the hour for the thermic and windshear development.

I’ll look at anything for those of us who do not know what we are doing, who lack the talent and are already so bad at this game they fly into hedges in flat fields a bit of a help up.

Leave a comment

XHTML: You can use these tags: <a href="" title=""> <abbr title=""> <acronym title=""> <blockquote cite=""> <code> <em> <strong>