Wednesday, 4 July 2018

postgis - Total great circle distance along a path of points


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.



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.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...