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.
I return to this issue because it is very similar to How do I find vector line bearing in QGIS or GRASS? and it can be solved with Python in the same way:
1) Haversine distance
One can find lots of scripts by searching Haversine distance with Python on the Internet and I choose one of them in Haversine Formula in Python (Bearing and Distance between two GPS points)
def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2, lat2])
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = math.sin(dlat/2)**2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon/2)**2
c = 2 * math.asin(math.sqrt(a))
km = 6367 * c
return km
We have a series of lines (points) in the file that must be treated in pairs (point1 - point2) to calculate the distance. For this we will use a simple iterator from Most pythonic way to get the previous element
def offset(iterable):
prev = None
for elem in iterable:
yield prev, elem
prev = elem
Now it is possible to read the file (example of Kerrie) in pairs of lines/points
import csv
with open('testhavers.csv', 'rb') as f:
reader = csv.DictReader(f)
for pair in offset(reader):
print pair
(None, {'LAT': '10.08527', 'LON': '124.59833', 'ID': '1', 'TIME': '21:24:37'})
({'LAT': '10.08527', 'LON': '124.59833', 'ID': '1', 'TIME': '21:24:37'},
{'LAT': '10.08523', 'LON': '124.59830', 'ID': '2', 'TIME': '21:25:07'})
({'LAT': '10.08523', 'LON': '124.59830', 'ID': '2', 'TIME': '21:25:07'},
{'LAT': '10.08526', 'LON': '124.59832', 'ID': '3', 'TIME': '21:25:37'})
({'LAT': '10.08526', 'LON': '124.59832', 'ID': '3', 'TIME': '21:25:37'},
{'LAT': '10.08526', 'LON': '124.59831', 'ID': '4', 'TIME': '21:26:07'})
Then create a shapefile containing the original fields of the csv file and a new field for the distance with the Python modules Shapely and Fiona of Sean Gillies:
import fiona
from shapely.geometry import Point, mapping
# creation of the schema of the shapefile (geometry and fields)
schema = { 'geometry': 'Point', 'properties':{'ID': 'int', 'LAT':'float', 'LON':'float', 'TIME':'str','distance' : 'float'}}
# creation of the shapefile:
with fiona.collection("result.shp", "w", "ESRI Shapefile", schema) as output:
# reading the csv file
with open('testhavers.csv', 'rb') as f:
reader = csv.DictReader(f)
# we need here to eliminate the first pair of point with None
for i, pair in enumerate(offset(reader)):
if i == 0: (pair with None)
# writing of the point geometry and the attributes
point = Point(float(pair[1]['LON']), float(pair[1]['LAT']))
dist = 0 # None
output.write({'properties': {'ID':int(pair[1]['ID']),'LAT':float(pair[1]['LAT']),'LON':float(pair[1]['LON']), 'TIME':pair[1]['TIME'],'distance': dist},'geometry': mapping(point)})
else:
# writing of the point geometry and the attributes
point = Point(float(pair[1]['LON']), float(pair[1]['LAT']))
# Haversine distance between pairs of points
dist = haversine(float(pair[0]['LON']), float(pair[0]['LAT']), float(pair[1]['LON']),float(pair[1]['LAT']))
output.write({'properties': {'ID':int(pair[1]['ID']),'LAT':float(pair[1]['LAT']),'LON':float(pair[1]['LON']), 'TIME':pair[1]['TIME'],'distance': dist},'geometry': mapping(point)})
and the result:
It is also possible to do it with PyQGIS but it is more complex than Fiona which uses simple dictionaries to create shapefiles.
You can use another function to calculate the Haversine distance (Why is law of cosines more preferable than haversine when calculating distance between two latitude-longitude points?) without any problem, only the distance calculation changes, not the process of creating the shapefile.
Best Answer
So your GPS coordinates are far apart? If they are close - eg. short segments along a road, then "straight line" is good enough.
There is an open source extension for PostGIS that can calculate road distances. It is called
pgRouting. You will also need to create a road database - typically OpenStreetMaps is used.