Scratching Surfaces

I returned from a short bike outing with my ride captured as a GPS track. Along the way, I also grabbed the rest stops as waypoints. Both of these were downloaded from the GPS as *.gpx files. So I have tracks.gpx and waypoints.gpx. Now I want to push these layers straight into Spatialite, and do some calculations.

I pulled the data from my GPS using gpsbabel but you can also choose the Quantum GIS GPSTools plugin to download GPS data in GPX format. Now to move these layers straight into a Spatialite database I use the ogr2ogr utility. This is part of the GDAL toolset, the “Swiss Army Knife” of spatial data formats. It allows me to move vector data from one format to another, with several options to filter or alter the original data along the way.
Here’s how I start:

$ ogr2ogr -f SQLite -nln ride -dsco "SPATIALITE=YES" -t_srs "EPSG:4326" bike_rides.sqlite bike_rides.gpx tracks

The first option sets the output format ‘-f‘ as SQLite. Next I use ‘-nln‘ to choose a “New Layer Name” for the table to be created. The ‘-dsco‘ options means “DataSet Create Option” and here we dictate Spatialite. The next flag ‘-t_srs‘ sets the spatial reference system for the data layer. While not strictly necessary (GPS data is always Lon/Lat WGS84) it’s a good policy to always set the SRS explicitly. Now pay attention to the order of the next parameters: the target database is first, and the source (gpx file) is second. This is backwards from many other command line tools. Normally we expect first the “from” followed by the “to”. For ogr2ogr it’s: “destination” first and “source” following. And finally at the end of the above command we choose to bring only GPS tracks into the database.

If our GPX file contains several tracks (perhaps from previous bike rides) we can filter those out with a simple -sql option such as:

$ ogr2ogr -f SQLite bike_rides.sqlite bike_rides.gpx -nln amazia_ride -dsco "SPATIALITE=YES" -sql "SELECT name, number FROM tracks WHERE name='Amazia-2' " tracks


Next, I need to get the stops along the way, which were captured as waypoints. Similar to above:

$ ogr2ogr -f SQLite bike_rides.sqlite waypoints.gpx -nln stops -append waypoints

Now, I don’t user the ‘-dsco’ parameter since the database already exists. What’s more, I must specify ‘-append‘ to instruct ogr2ogr to add this layer as a new table in the existing sqlite database.

OK, how long was the whole ride? The built in Spatialite function GeodesicLength() calculates lengths of Lon/Lat features in meters, along the WGS84 ellipsoid.

spatialite> SELECT GeodesicLength(r.Geometry) AS “Total Length” FROM ride AS r;
Total Length
—————
9696.2803124924

But I want to know how much each section of the ride was. In other words, I want to see the accumulated distance at each of the stops. We have in spatialite a function ST_Line_Locate_Point() which returns the distance along a line feature for points near the line. (It actually finds the nearest point on the line for each of the input points). So we do:

spatialite> SELECT s.name AS "Stop name",
...> ST_Line_Locate_Point(r.Geometry, s.Geometry) AS "Stop Locations"
...> FROM ride AS r, stops AS s
...> GROUP BY s.name ORDER BY s.OGC_FID;
Stop name   Stop Locations
----------  --------------
start       0.0
rest 1      0.124783450768
water stop  0.361703783084
rest 2      0.562004941906
rest 3      0.712163267248
overlook    0.850818111308
end         0.996608235043


The OGC_FID field is automatically created by ogr2ogr as a primary key. So doing “ORDER BY” on that field insures that the stop points are displayed in the correct order. But the numbers above are not distances, but rather the fraction of the total line length where each point falls. Now multiplying these fractions by the total length from the first query above, I get:

spatialite> SELECT s.name AS "Stop name",
...> ST_Line_Locate_Point(r.Geometry, s.Geometry)*GeodesicLength(r.Geometry) AS "Accumlulated Distances"
...> FROM ride AS r, stops AS s
...> GROUP BY s.name ORDER BY s.OGC_FID;
Stop name   Accumlulated Distances
----------  ----------------------
start       0.0
rest 1      1209.93531700926
water stop  3507.18127087581
rest 2      5449.35745373272
rest 3      6905.3346675029
overlook    8249.77090219153
end         9663.3928087156

Nice. One last step, let’s put these distances permanently into the “stops” table.

spatialite> ALTER TABLE stops ADD COLUMN accum_dist float;
spatialite> UPDATE stops SET accum_dist=(SELECT
...> ST_Line_Locate_Point(r.Geometry, s.Geometry)*GeodesicLength(r.Geometry)
...> FROM tracks AS r, stops AS s
...> WHERE s.OGC_FID=stops.OGC_FID);

We use a subquery to UPDATE the stops table. Take note of the WHERE condition: the OGC_FID values from the stops table outside of the subquery must match the values from the same column inside the query.

spatialite> SELECT name AS "Stop name",accum_dist AS "Accumulated Distance" FROM stops;
Stop name Accumulated Distance
---------- --------------------
start       0.0
rest 1      1209.93531700926
water stop  3507.18127087581
rest 2      5449.35745373272
rest 3      6905.3346675029
overlook    8249.77090219153
end         9663.3928087156

There we are: accumulated distances for waypoints along a GPS track.

This work, unless otherwise expressly stated, is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.