Freesteel Blog

GPS is a jerk

Monday, March 20th, 2017 at 1:45 pm Written by: Julian

Last week I finally had my first flight of the year with my newly build flight data logger. I can’t believe the number of issues it’s already thrown up.

At least I may be making quick enough progress to get past the issues (rather than being swamped by them) using this exceptionally efficient Jupyter/Pandas technology.

For example, my code for parsing and loading the IGC file is 15 lines long.

The code for loading in my flight logger data into a timeseries is just as brief, if you consider each data type individually (there are more than 13 of them from humidity sensors to an orientation meter).

The GPS time series from my flight logger (at 10Hz) can be crudely converted it to XYs in metres, like so:

# pQ is the GPS position pandas.DataFrame
earthrad = 6378137
lng0, lat0 = pQ.iloc[0].lng, pq.iloc[0].lat
nyfac = 2*math.pi*earthrad/360
exfac = nyfac*math.cos(math.radians(lat0))
pQ["x"] = (pQ.lng - lng0)*exfac
pQ["y"] = (pQ.lat - lat0)*nyfac
plt.plot(pQ.x, pQ.y)

Note the suspicious sharp turn near (-1000, -400). Here’s another sharp turn somewhere else in the sequence covering a 1 minute 5 second period using time slicing technology:

t0, t1 = Timestamp("2017-03-09 15:42:55"), Timestamp("2017-03-09 15:44:00")
q = fd.pQ[t0:t1]
plt.plot(q.x, q.y)

The dot is at the start point, time=t0.

It’s taken me days to conclude that this is complete bollocks.

The first piece of evidence was the summation of the velocities made by the GPVTG records of velicity and heading (degrees), which we can sum cumulatively like this to recreate the path:

v = fd.pV[t0:t1]
vsx = (v.vel*numpy.sin(numpy.radians(v.deg))).cumsum()
vsy = (v.vel*numpy.cos(numpy.radians(v.deg))).cumsum()
plt.plot(vsx*0.1, vsy*0.1)  # *0.1 because the readings are 10Hz

No impossibly sharp U-turns there, but it does claim I did an S-turn, first to the left to head south before turning sharp right to turn north. At no time does this show me heading east.

This was not consistent with the BNO055 orientation sensor, which reads at 100Hz and which can be summed in the same way as the GPS velocity on the assumption that the glider heads at 10m/s in the direction it is pointed. (Don’t worry about the wind. It was very light that day.)

dorient = utils.NorthOrient(fd.pZ)[t0:t1]
vosx = numpy.sin(numpy.radians(dorient)).cumsum()
vosy = numpy.cos(numpy.radians(dorient)).cumsum()
plt.plot(vosx*0.1, vosy*0.1)

This is a lot more realistic, as well as exhibiting a 270 degree left hand turn not shown by the GPS.

It is, however, consistent with the video footage from my badly misaligned keel camera.

The GPS does not always produce bollocks. How do we find out when it is good?

We can calculate the heading in degrees from consecutive locations and compare it to the velocity heading to observe that sometimes they do agree and sometimes they are very far out, like so:

q = fd.pQ[t0:t1]
vq = q.shift(5) - q.shift(-5)  # 1-second range at 10Hz
qdeg = (numpy.degrees(numpy.arctan2(vq.x, vq.y))+180)
qdeg.plot()
fd.pV[t0:t1].deg.plot()

The same can be done in the velocity domain:

numpy.sqrt(vq.x**2 + vq.y**2).plot()
v.vel.plot()

And here’s how it resolves itself in the vy domain, which is going to be easier to compare:

vx = v.vel*numpy.sin(numpy.radians(v.deg))
vy = v.vel*numpy.cos(numpy.radians(v.deg))
vq = q.shift(-5) - q.shift(5)
vy.plot()
vq.y.plot()

The readings and derived readings are a bit hard to compare like this, so I’ve interpolated the velocities to the gps position time series (because they have slightly different time-stamps) and low-pass filtered them to get rid of some of the noise and width alignment that would interfere with the comparison.

vyi = utils.InterpT(vq.y, vy)
plt.plot(utils.FiltFiltButter(vyi, 0.015), label="posvel")
plt.plot(utils.FiltFiltButter(vq.y, 0.015), label="realvel")

What we would like to work up to is a method of selecting out the sets of points in the gps sequence that are bad so we don’t waste time attempting to process bad data.

Here’s how we plot the sequences of points (in black) where the gps position and velocity disagree by more than 4m/s in either dimension.

# from the position, get the velocity in x and y
q = fd.pQ[tf0:tf1]  
vq = (q.shift(-5) - q.shift(5)).dropna()
    # drop the 5 end-values as they interfere with the filter

# from the velocity calculate the x and y
v = fd.pV[tf0:tf1]
vx = v.vel*numpy.sin(numpy.radians(v.deg))
vy = v.vel*numpy.cos(numpy.radians(v.deg))

# interpolate to position series and filter the difference
vxi = utils.InterpT(vq.x, vx)
vxd = numpy.abs(utils.FiltFiltButter(vxi - vq.x, 0.015))
vyi = utils.InterpT(vq.y, vy)
vyd = numpy.abs(utils.FiltFiltButter(vyi - vq.y, 0.015))

# *advanced pandas code here*
# drop the 5 end-values to allow comparison
# and subselect against boolean timeseries specified by when the 
# difference in smoothed velocities exceeds 4 in x or y axis
badq = q.iloc[5:-5][(vxd>4) | (vyd>4)]
plt.plot(q.x, q.y)
plt.scatter(badq.x, badq.y)

And this was when I remembered that I was carrying a second GPS in my 6030 vario device and had to quickly write a parser for its IGC file, to plot it as an overlay like so:

So it’s plausible that by comparing the GPS position differences to the GPS velocity we’re able to spot the bits where the GPS readings are pants so as not to waste any effort trying to make sense of it, because there is no sense.

For a closer look, here’s how they line up in that 1 minute period shown above:

How about comparing these two GPS readings to each other?

The vario GPS reads once every 2 seconds, while my flight logger GPS reads 10x a second, but with a bit of pandas code, we can align GPS timestamps u (rather than by the microcontroller’s millisecond timestamp) and subselect accordingly:

useries = g.index.to_series().apply(lambda X: X.hour*3600+X.minute*60+X.second+0.0)
dfg = pandas.DataFrame(data={"lat":g.lat, "lng":g.lng, "u":useries})
dfg = dfg.set_index("u")

dfg["qlat"] = pandas.Series(list(q.lat), q.u/1000)
dfg["qlng"] = pandas.Series(list(q.lng), q.u/1000)

plt.scatter((dfg.lat - dfg.qlat)*50, (dfg.lng - dfg.qlng)*50)

It’s a donut!

There’s something suspicious about donut errors because it looks like something is chasing its tail.

Let’s try adding on 2 seconds to the datalogger data to see how it compares:

qu = q.u/1000+2.0
dfg["qx"] = pandas.Series(list(q.x), qu)
dfg["qy"] = pandas.Series(list(q.y), qu)
plt.scatter((dfg.x - dfg.qx), (dfg.y - dfg.qy))

And if we add on 4.0 seconds we get the donut back:

We can loop through to find the offset which has the minimal error, like so:

varss = [ ]
for i in range(40):
    qu = q.u/1000+i/10
    dfg["qx"] = pandas.Series(list(q.x), qu)
    dfg["qy"] = pandas.Series(list(q.y), qu)
    varss.append((dfg.x - dfg.qx).var() + (dfg.y - dfg.qy).var())
plt.plot(varss)

Turns out this is at 2.1seconds.

What basically happens is that the timestamps in the vario match the position given in the timestamps by the flight logger two seconds earlier.

I don’t know if this (or the horrendous directional errors) are caused by setting the GPS to sample at such a high rate (10 times a second), or it being a comparatively bad device, or what.

It appears I can’t trust the GPS readings to even log my direction to within 90degrees accuracy, so there’s no point in using it in relation to the glider kinematics (turning circles etc). This explains why all of my orientation kinematics attempts have fallen flat.

The GPS velocity method for sifting out bollocks sets of readings is plausible, but probably too noisy and messy.

Given that different GPS devices can at times disagree, it looks like I could bolt on one or two extra GPS chips into my unit and maybe run them at a different update frequencies. And then only run kinematics models during periods of time when they all agree.

To recap, the challenge is to discover a device that can inform you, the pilot, of three things in real-time: (1) the glider’s performance, (2) your own performance, and (3) important facts about the airmass. This should make a difference for people like me who don’t have the talent and good fortune to be able to practice flying for hundreds of hours per year.

Next up, air temperature, sunlight and very shoddy barometer readings.

Update

A long afternoon of bodging code, resoldering wires to access serial pins and whole stick of hot glue and I now have a second GPS stuck on to the top of the system.

Nobody has ever accused me of wasting too much time on design, but you’ve got to keep focused on the problem.

There is still scope to add on a third GPS chip should it be required, but this will hopefully not be necessary. The point is not to have an accurate position all the time, it’s to know when the position is not at all accurate so I don’t attempt to base any calculations on it when it is bad.