# How to interpolate GPS Positions in PostGIS

I have a PostGIS table of GPS positions for every five seconds:

``````2011-01-01 00:00:05, POINT(x1,y1)
2011-01-01 00:00:10, POINT(x2,y2)
2011-01-01 00:00:15, POINT(x3,y3)
...
``````

I'm looking for a query that will return values (timestamp and point) for every second. It's ok to assume that points are connected by a straight line.

I'm specifically looking for a way to do this inside the database and not by writing some external script.

-
I think that you will need to write a PL/Python function for that. – Pablo May 16 '11 at 19:23
Here is a snippet from postgis in action that may help: bostongis.com/postgis_translate.snippet – Pablo May 16 '11 at 19:34
@Pablo: Yes, most likely. I'll adjust my question. – underdark May 16 '11 at 19:35

hallo

If your original table is called gps_p, your timestamp field is called ts and the points is called th_geom:

``````SELECT (geom).geom,  ts1 + (((geom).path[1]-1) ||' seconds')::interval FROM
(SELECT ts1, ST_DumpPoints(ST_Segmentize(geom, ST_Length(geom)/5)) as geom FROM
(SELECT ts1, ST_LineFromMultipoint(ST_Union(geom1, geom2)) as geom FROM
(SELECT p1.ts as ts1, p2.ts as ts2, p1.the_geom as geom1, p2.the_geom as geom2
FROM gps_p p1 INNER JOIN gps_p p2 on p1.ts + '00:00:05'::interval = p2.ts
) a
)b
) c
WHERE (geom).path[1] <= 5;
``````

What it does is that it builds lines between the points and use st_segmentize to divide the line in 5 segments.

If it is not exactly 5 seconds between your original points it will not work. Then you can just add an id field with a sequence and use that to selfjoin the table with id1+1 = id2 instead.

HTH

/Nicklas

-

here is a code draft for pl/python, it's just the basic idea of translating the points by a given distance and azimuth.
To run postgis functions in pl/python the only solution that I found is to use plpy.prepare and plpy.execute (very boring).

``````total_distance=St_distance(P1,P2)
azimuth=st_azimuth(p1,p2)
partial_distance=total_distance / 5

for i in range(4):
distance = (i+1)*partial_distance
x_increment=distance*math.cos(math.degrees(azimuth))
y_increment=distance*math.sin(math.degrees(azimuth))
ST_translate(P1, x_increment, y_increment)
``````
-

If I am not mistaken...
What you would need to do is to determine the connecting line and then do a divide on it.

-