gps

So, on the bike, I use a portable GPS to keep track of my speed and to track the mileage done on the bike so I know when to next put it in for service. Originally I just relied on the trip counter in the GPS, but then found that this could develop quite an error if left to tick over for a few months.

Thus, I wrote a simple CGI application in Perl and SQLite3 that would track the odometer readings. Plain, simple, and it’s worked quite well, but remembering to punch in the current odometer reading is a chore, and my stats are only as granular as I submit: if I want to see what distance I did on a particular day, I either have to have had the foresight to store readings at the start and end of that day, or I’m stuffed.

I also keep the GPX tracklogs. While the Garmin 650 is not great at handling lots of tracklogs (and for some moronic reason, they name the files “Day DD-MMM-YY HH.MM.SS.gpx”, not something sensible like “YYYY-MM-DDTHH-MM-SS.gpx”), it’s good enough that I can periodically siphon off the track logs for storage on my laptop. I then have a record of where I’ve been.

Theoretically, this also has the distance travelled, I could make a service that just consumes the GPX files, and tallies up the distances that way. Maybe even visualise heat-map style, where I go most. (No prizes for guessing “work” … but where else?)

The existing system uses SQLite, and specifically, its views, as poor man’s stored procedures. It’s hacky, inefficient, and sooner or later I’ll have performance problems. PostGIS is an extension onto PostgreSQL which supports a large number of spatial operations, including finding the length of a series of points, which is exactly the problem I’m trying to solve right now. The catch is, how do you import the data?

Enter GDAL

GDAL is a library of geographic functions for answering these kinds of questions. It ships with a utility ogr2ogr, which can take geographic information in a variety of formats, and convert to a variety of output formats. Crucially, this tool supports consuming GPX files and writing to a PostGIS database.

Loading one file, is easy enough:

$ ogr2ogr -oo GPX_ELE_AS_25D=YES \
  -dim 3 \
  -gt 65536 \
  -lco GEOM_TYPE=geography \
  -preserve_fid \
  -f PostgreSQL \
  "PG:dbname=yourdb" yourfile.gpx \
  tracks track_points

The arguments here were found by trial-and-error.  Specifically, -oo GPX_ELE_AS_25D=YES and -dim 3 tell ogr2ogr to preserve the elevation in the point information (as well as keeping a copy of it in the ele column). -lco GEOM_TYPE=geography tells ogr2ogr to use the geography data type in PostGIS.

Look in the database, and you’ll see two tables, tracks and track_points. Sadly, you don’t get to choose the names of these (not easily anyway, there is -nln, but it then will create one table with the given name, put the tracks in it, then blow it away and replace it with a table of the same name containing points), and there’s no foreign keys between the two.

The fun starts when you try to import a second GPX file. Run that command again, and because of -preserve_fid, you’ll get a primary key clash. Take that away, and the track_fid column in track_points becomes meaningless.

If you drop -preserve_fid, then track_fid gets set to 0 for all points.  Useless.

Importing many GPX files

Out of the box, this just wasn’t going to fly, so we needed to do things a little different.  Firstly, I duplicated the schema that GDAL creates, creating my own tables which will ultimately store the data.  I then used a wrapper shell script that calls psql before and after ogr2ogr so I can re-map the primary keys to maintain relationships.

Schema SQL

CREATE SEQUENCE public.gpx_points_ogc_fid_seq
    INCREMENT 1
    START 1
    MINVALUE 1
    MAXVALUE 2147483647
    CACHE 1;

CREATE SEQUENCE public.gpx_tracks_ogc_fid_seq
    INCREMENT 1
    START 1
    MINVALUE 1
    MAXVALUE 2147483647
    CACHE 1;

CREATE TABLE public.gpx_tracks
(
    ogc_fid integer NOT NULL,
    name character varying COLLATE pg_catalog."default",
    cmt character varying COLLATE pg_catalog."default",
    "desc" character varying COLLATE pg_catalog."default",
    src character varying COLLATE pg_catalog."default",
    link1_href character varying COLLATE pg_catalog."default",
    link1_text character varying COLLATE pg_catalog."default",
    link1_type character varying COLLATE pg_catalog."default",
    link2_href character varying COLLATE pg_catalog."default",
    link2_text character varying COLLATE pg_catalog."default",
    link2_type character varying COLLATE pg_catalog."default",
    "number" integer,
    type character varying COLLATE pg_catalog."default",
    gpxx_trackextension character varying COLLATE pg_catalog."default",
    the_geog geography(MultiLineStringZ,4326),
    CONSTRAINT gpx_tracks_pkey PRIMARY KEY (ogc_fid)
)
WITH (
    OIDS = FALSE
)
TABLESPACE pg_default;

CREATE TABLE public.gpx_points
(
    ogc_fid integer NOT NULL,
    track_fid integer,
    track_seg_id integer,
    track_seg_point_id integer,
    ele double precision,
    "time" timestamp with time zone,
    magvar double precision,
    geoidheight double precision,
    name character varying COLLATE pg_catalog."default",
    cmt character varying COLLATE pg_catalog."default",
    "desc" character varying COLLATE pg_catalog."default",
    src character varying COLLATE pg_catalog."default",
    link1_href character varying COLLATE pg_catalog."default",
    link1_text character varying COLLATE pg_catalog."default",
    link1_type character varying COLLATE pg_catalog."default",
    link2_href character varying COLLATE pg_catalog."default",
    link2_text character varying COLLATE pg_catalog."default",
    link2_type character varying COLLATE pg_catalog."default",
    sym character varying COLLATE pg_catalog."default",
    type character varying COLLATE pg_catalog."default",
    fix character varying COLLATE pg_catalog."default",
    sat integer,
    hdop double precision,
    vdop double precision,
    pdop double precision,
    ageofdgpsdata double precision,
    dgpsid integer,
    the_geog geography(PointZ,4326),
    CONSTRAINT gpx_points_pkey PRIMARY KEY (ogc_fid),
    CONSTRAINT gpx_points_track_fid_fkey FOREIGN KEY (track_fid)
        REFERENCES public.gpx_tracks (ogc_fid) MATCH SIMPLE
        ON UPDATE RESTRICT
        ON DELETE RESTRICT
)
WITH (
    OIDS = FALSE
)
TABLESPACE pg_default;

The wrapper script

 1 #!/bin/sh
 2 
 3 DB=tracklog
 4 
 5 for f in "$@"; do
 6         psql tracklog <<EOF
 7 DROP TABLE IF EXISTS tracks;
 8 DROP TABLE IF EXISTS track_points;
 9 EOF
10         ogr2ogr -oo GPX_ELE_AS_25D=YES \
11                 -dim 3 \
12                 -gt 65536 \
13                 -lco SPATIAL_INDEX=FALSE \
14                 -lco GEOM_TYPE=geography \
15                 -overwrite \
16                 -preserve_fid \
17                 -f PostgreSQL \
18                 "PG:dbname=${DB}" "$f" \
19                 tracks track_points
20 
21         # Re-map FIDs then insert into real tables.
22         psql tracklog <<EOF
23         CREATE TEMPORARY TABLE track_fids AS
24         SELECT  ogc_fid AS orig_fid,
25                 nextval('gpx_tracks_ogc_fid_seq') AS ogc_fid
26         FROM    tracks;
27 
28         CREATE TEMPORARY TABLE point_fids AS
29         SELECT  ogc_fid AS orig_fid,
30                 nextval('gpx_points_ogc_fid_seq') AS ogc_fid
31         FROM    track_points;
32 
33         INSERT INTO gpx_tracks
34         SELECT  track_fids.ogc_fid AS ogc_fid,
35                 tracks.name as name,
36                 tracks.cmt as cmt,
37                 tracks."desc" as "desc",
38                 tracks.src as src,
39                 tracks.link1_href as link1_href,
40                 tracks.link1_text as link1_text,
41                 tracks.link1_type as link1_type,
42                 tracks.link2_href as link2_href,
43                 tracks.link2_text as link2_text,
44                 tracks.link2_type as link2_type,
45                 tracks."number" as "number",
46                 tracks.type as type,
47                 tracks.gpxx_trackextension as gpxx_trackextension,
48                 tracks.the_geog as the_geog
49         FROM    track_fids, tracks
50         WHERE   track_fids.orig_fid=tracks.ogc_fid;
51 
52         INSERT INTO gpx_points
53         SELECT  point_fids.ogc_fid AS ogc_fid,
54                 track_fids.ogc_fid AS track_fid,
55                 track_points.track_seg_id AS track_seg_id,
56                 track_points.track_seg_point_id AS track_seg_point_id,
57                 track_points.ele AS ele,
58                 track_points."time" AS "time",
59                 track_points.magvar AS magvar,
60                 track_points.geoidheight AS geoidheight,
61                 track_points.name AS name,
62                 track_points.cmt AS cmt,
63                 track_points."desc" AS "desc",
64                 track_points.src AS src,
65                 track_points.link1_href AS link1_href,
66                 track_points.link1_text AS link1_text,
67                 track_points.link1_type AS link1_type,
68                 track_points.link2_href AS link2_href,
69                 track_points.link2_text AS link2_text,
70                 track_points.link2_type AS link2_type,
71                 track_points.sym AS sym,
72                 track_points.type AS type,
73                 track_points.fix AS fix,
74                 track_points.sat AS sat,
75                 track_points.hdop AS hdop,
76                 track_points.vdop AS vdop,
77                 track_points.pdop AS pdop,
78                 track_points.ageofdgpsdata AS ageofdgpsdata,
79                 track_points.dgpsid AS dgpsid,
80                 track_points.the_geog AS the_geog
81         FROM    track_points, track_fids, point_fids
82         WHERE   point_fids.orig_fid=track_points.ogc_fid
83         AND     track_fids.orig_fid=track_points.track_fid;
84 
85         DROP TABLE tracks;
86         DROP TABLE track_points;
87         DROP TABLE track_fids;
88         DROP TABLE point_fids;
89 EOF
90 done

Getting the length of a track

Having imported all the data, we can do something like this:

SELECT ogc_fid, name,
  ST_Length(the_geog, false)/1000 as dist_in_km
FROM gpx_tracks order by ogc_fid desc limit 10;

and get this:

Visualisation with QGIS

Turns out, this is straightforward…

1. In your workspace, there’s a tree with the different layer types you can add, including PostGIS… right-click on this and select New Connection… fill in the details for your PostgreSQL database.
2. Below that is XYZ Tiles…, right click again, select New Connection for OpenStreetMap, and use the URL https://a.tile.openstreetmap.org/{z}/{x}/{y}.png (also, see their policy).
3. Drag the OpenStreetMap connection to your layers
4. Expand the PostGIS connection you just made, and look for the gpx_tracks table, drag this on top of your OpenStreetMap layer.

Below is everywhere I’ve been with the GPS tracklog running.  Much of what you see is the big loop a few of us did in 2012, including my trip to Ballarat for the 2012 LCA.

If I zoom in on Brisbane, unsurprisingly, some areas show up very clearly as being common haunts for me:

A bit of SQL voodoo, and I come up with this:

In orange is the territory covered on the Boulder (minus what was covered before I got the GPS), in blue the territory covered on the Talon 29 ER 0, and in red, on my current commuter (Toughroad SLR2).