[GIS] Total great circle distance along a path of points

distancepointpostgis

I am looking for an SQL query that sums the distances between a series of points using PostGIS.

The situation is as follows:

I have 2 tables: tracks and points.
Points belong to a track and thus each points record has a track_id field.
Points also have geom field (in WGS84 = GPS coordinates) and a timestamp.

What I would like to do is calculate the total (geographically correct) distance for each track.

In orther words, I need a query that, when given a certain track_id, looks up all the points for that track, orders them by timestamp and than calculates the distance between each consecutive pair of points and sums these distances to give me the total distances of the track in question.

Best Answer

It is fairly straightforward to write a single query for this, you just need to break it down into logical blocks. To start with, you need your point geometry sorted by track and time:

SELECT track_id, gpstime, ST_AsText(the_geom) FROM points ORDER BY track_id, gpstime;

Here I assume you have your points in a geometry column. A geography column in this case would achieve the same thing, and actually save a bit of faff in the next step.

(Note here I've used ST_AsText() to make things readable. You should remove this in the final query).

With my simple database of 2 tracks and 6 points, that query gives me:

 track_id |            gpstime            |         st_astext
----------+-------------------------------+----------------------------
        1 | 2012-01-27 14:36:26.354665+00 | POINT(-71.064544 42.28787)
        1 | 2012-01-27 14:36:53.608661+00 | POINT(-71.063544 42.28987)
        1 | 2012-01-27 14:37:10.841856+00 | POINT(-71.066544 42.3)
        2 | 2012-01-27 14:37:46.400954+00 | POINT(-0.8 53.2)
        2 | 2012-01-27 14:37:55.752122+00 | POINT(-0.81 53.21)
        2 | 2012-01-27 14:38:07.773077+00 | POINT(-0.82 53.22)

The next step is to aggregate those points into linestrings. This is done with the ST_MakeLine() function:

SELECT gps.track_id, ST_AsText(ST_MakeLine(gps.the_geom)) AS track_line
    FROM (SELECT track_id, gpstime, the_geom FROM points ORDER BY track_id, gpstime) as gps
    GROUP BY gps.track_id
    ORDER BY gps.track_id;

Which yields:

 track_id |                             track_line
----------+---------------------------------------------------------------------
        1 | LINESTRING(-71.064544 42.28787,-71.063544 42.28987,-71.066544 42.3)
        2 | LINESTRING(-0.8 53.2,-0.81 53.21,-0.82 53.22)

Finally you need to get the length of that line over the spheroid of your choice. If you're using GPS then inevitably it'll be the WGS84 spheroid. The function to use is ST_Length_Spheroid().

So putting that all together, my final query ends up looking like this:

SELECT gps.track_id, tracks.name, ST_Length_Spheroid(ST_MakeLine(gps.the_geom), 'SPHEROID["WGS 84",6378137,298.257223563]') AS track_len
    FROM (SELECT track_id, gpstime, the_geom FROM points ORDER BY track_id, gpstime) as gps, tracks
    WHERE gps.track_id = tracks.track_id
    GROUP BY gps.track_id, tracks.name
    ORDER BY gps.track_id;

And I get these results:

 track_id | name  |    track_len
----------+-------+------------------
        1 | Larry | 1389.08015675763
        2 | Curly | 2596.09131911171

If you want to limit it to a specific track, change your WHERE clause to something like this:

 WHERE gps.track_id = tracks.track_id AND tracks.name = 'Larry'

If you're using geography columns instead of geometry, you can eliminate the spheroid part and just use ST_Length(). Note that the second parameter of ST_Length() when using geography columns should be set to TRUE if you want more accurate but slower calculations. The only drawback is that geography columns are limited to the WGS84 spheroid.

Related Question