Separating disjoint GPS path segments in Postgres

Question

This involves a lot of code; if it's better suited for StackOverflow, please let me know!

I'm working on a project involving trips that are constructed from GPS coordinates. I'm mapping the recorded trips on a webpage, with version 0.35 of Windshaft. Trips are collected from a mobile application, and the GPS data is inserted into a Postgres (9.3) database (with PostGIS (2.1.2) installed). The relevant tables are:

coord_geog:

 Column  |            Type             |
---------+-----------------------------|
id       | integer                     | 
trip_id  | integer                     |
recorded | timestamp without time zone |
geog     | geography(Point,4326)       |
geom     | geometry(Point,4326)        |
next     | integer                     |

trip_geom:

 Column  |            Type             |
---------+-----------------------------|
id       | integer                     |
purpose  | character varying(255)      |
start    | timestamp without time zone |
stop     | timestamp without time zone |
geom     | geometry(LineString,4326)   |

The trip_id column of coord_geog corresponds to the id column of trip_geom. The application supplies data for the id, trip_id, recorded, and geog columns of coord_geom, and the id, purpose, start, and stop columns of trip_geom. I'm actually not sure what the geom column of coord_geog is used for, but I'm too scared to remove it (and it doesn't appear to be used in anything mentioned in this question).

I'm not very sure how the whole Windshaft thing works, but if it helps, the Windshaft configuration is currently set as:

req.params.sql = "(select * from trip_geom_frag where purpose ilike '" + req.params.purpose + "%') as trip_geom_frag";
req.params =  _.extend({}, req.params, {style: style});

I'll specify trip_geom_frag lower down - it is similar to trip_geom.

The specific task I'm trying to improve is constructing a LineString path from all of the points in a path.

I inherited a Python script that's supposed to handle this task. Although I'm not particularly well-versed in Python or SQL, I recognized that a lot of things that the script was doing could be done as pure SQL, which resulted in a major speedup at first. After the points are anonymized (trimmed from each end of a trip), a query like this constructs the LineString:

UPDATE trip_geom t SET geom=(SELECT ST_MakeLine(line.geom) FROM (SELECT c.recorded, c.geom FROM coord_geog c WHERE c.trip_id=t.id ORDER BY c.recorded ASC) as line) WHERE t.id=%s;

where the %s parameter is an id field from a row in trip_geom. This constructs a LineString connecting the points in the trip. However, there tend to be errors in the GPS data, where a bad GPS fix causes one or more points to be incorrectly located far away, or the recording is paused and later resumed, resulting in a large 'jump' in the trip.

Valid points are generally close together (since the trips are bicycle trips at a relatively low speed), so when a trip has all valid GPS data, the LineString is a fairly smooth plot (map background by CartoDB):

However, many trips (especially the long ones) have jumps or errors, resulting in trips that look like this:

While I'm sure there are some highly accomplished bicyclists, I doubt that someone was able to bike straight across a river; the straight line appears to be caused by the bicyclist pausing and resuming a trip recording. These errors aren't a big deal when isolated, but when there are many of them, the map of trips becomes very cluttered with straight lines that cross over large areas of the map.

Therefore, the specific task I am attempting to accomplish is to generate LineString paths without adjacent points that are very far apart. I imagine that this is a task well suited to PostGIS (or SQL in general), but I'm having difficulty finding a solution.

My first attempted solution was to remove points that were far from the previous point. Since I don't have the exact SQL any more, here is the procedure in pseudo-code:

for t in trips:
    i = 1
    while i < t.points.length:
        current_pt = t.points[i]
        last_pt = t.points[i-1]
        if ! ST_DWithin( current_pt, last_pt, 100 ):
            delete_point( current_pt )
        else:
            i = i + 1

This resulted in smooth trips, but the trips with jumps in the GPS data were truncated after the jump. For example, in a trip with a jump between the trip's first and second points, all of the points in the trip after the first point would be removed.

My second (and current) solution is to create a second table like trip_geom to hold LineStrings from trips, which are separated whenever there's a jump between points. I decided to call these fragmented paths - there's probably a better term. Here's the new table I created:

trip_geom_frag:

  Column |           Type              |
---------+-----------------------------|
id       | integer                     |
geom     | geometry(LineString,4326)   |
purpose  | character varying(255)      |
orig_trip| integer                     |

Then, using the next column of coord_geog (which I added), I use a query like this to construct a sort of linked list of the points in each trip, so that a -1 value represents either a jump in the points or the end of a trip:

UPDATE coord_geog cur_pt SET next = CASE WHEN ST_DWithin( cur_pt.geog, (SELECT geog FROM coord_geog WHERE trip_id=cur_pt.trip_id AND recorded > cur_pt.recorded ORDER BY recorded ASC LIMIT 1), %s) THEN 1 ELSE -1 END;

Then, in Python, I insert a row into trip_geom_frag for each segment of a trip, with the end of a segment being marked by a -1:

for t in trips:
    trip_id = t.id

    c.execute('SELECT id, next FROM coord_geog WHERE trip_id=%s ORDER BY id ASC;', (trip_id,));
    coords = c.fetchall()
    ind = [coord[1] for coord in coords]

    while len(coords) > 0:
        i = ind.index(-1)
        c.execute('INSERT INTO trip_geom_frag (geom, orig_trip, purpose) SELECT ST_MakeLine(line.geom), %s, %s FROM (SELECT c.geom FROM coord_geog c WHERE c.trip_id=%s AND c.id >= %s AND c.id <= %s ORDER BY c.recorded ASC) as line;', (trip_id, purpose, trip_id, coords[0][0], coords[i][0], ))
        del coords[:i+1]
        del ind[:i+1]

Again, I'm not very experienced in Python or SQL, so forgive me if this is horrible; I'm pretty sure that there is some way to do this without shuffling the data between Python and Postgres.

In the best case of no GPS jumps or errors, this produces a single LineString, or, in the case of jumps or errors, possibly many LineStrings. This produces paths that are separated whenever there's a jump between points (same data pictured):

This data is exactly like I'd hope it would come out, so this is a sufficient solution for the time being. However, there are currently 800 trips and approximately 1.5 million total points; processing them took about four days (on a crappy server) - I don't know if that's because there's something wrong with Postgres or if it's just because my Python solution is so inefficient. Either way, if for some reason I need to re-process all of the trips, it could take far too long before the trips are ready.

As I see it, here are a few approaches that may work well:

• Insert coordinates into the LineString in a way that allows particular segments to be excluded
• Partition each trip into contiguous segments, by changing the arguments to ST_MakeLine
• Make Windshaft ignore segments of a LineString that exceed some length
• Scrap ST_MakeLine, and instead make a custom path that somehow skips over a jump in points
• Something obvious that I don't know about because I'm too unfamiliar with the tools present

My hope is that there's some simple solution, ideally requiring only PostGIS, that can figure out where the points in a trip should be partitioned and insert them into a table. My question, therefore, is what approach should I take to partition paths that are used in calls to ST_MakeLine to only draw contiguous points?

Since I'm already saving the data as it was previously constructed (non-fragmented) in a separate table, I'm free to modify any of the tables or mangle any of the data if it helps in some way.

---

Update: I set up Postgres on my (much more powerful) local machine and downloaded the database, re-running the processing locally. It's much faster - it looks like all of the points can be processed in under 30 minutes, meaning that the culprit was indeed a crappy server. It would still be great to have a solution in pure SQL, though.

Answer 1

This is my attempt for a query that should create parts of trips, splitted at jumps of about 100 m.

CREATE TABLE trip_geom_parts (
    trip_id integer NOT NULL,
    part_id integer NOT NULL,
    start timestamp without time zone,
    stop timestamp without time zone,
    geom geometry(linestring, 4326),
    PRIMARY KEY (trip_id, part_id));
CREATE TEMPORARY SEQUENCE seq;
INSERT INTO trip_geom_parts (trip_id, part_id, start, stop, geom)
    WITH s1 AS (SELECT
        trip_id,
        CASE WHEN trip_id = lag(trip_id) OVER w
             THEN CASE WHEN ST_Distance_Sphere(geom, lag(geom) OVER w) < 100.0
                       THEN currval('seq')
                       ELSE nextval('seq')
                  END
             ELSE setval('seq', 1)
        END AS part_id,
        recorded,
        geom
        FROM coord_geog
        WINDOW w AS (PARTITION BY trip_id ORDER BY recorded))
    SELECT
        trip_id, part_id,
        min(recorded) AS start, max(recorded) AS stop,
        ST_MakeLine(geom ORDER BY recorded)
        FROM s1
        GROUP BY trip_id, part_id;

For better performance it may help to

CREATE INDEX ON coord_geog (trip_id, recorded);

I didn't test it, so may need some debugging.

Answer 2

Not directly related to your question, but sounds like you could find the Trackintel library useful: https://github.com/mie-lab/trackintel (Disclaimer: I'm part of the lab developing it)

It offers easy-to-use methods to load data from PostGIS and do many preprocessing steps such as extracting triplegs and trips, as well as further analysis (clustering etc), and visualisation.