So interpolieren Sie GPS-Positionen in PostGIS

13

Ich habe alle fünf Sekunden eine PostGIS-Tabelle mit GPS-Positionen:

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)
...

Ich suche nach einer Abfrage, die Werte (Zeitstempel und Punkt) für jede Sekunde zurückgibt. Es ist in Ordnung anzunehmen, dass Punkte durch eine gerade Linie verbunden sind.

Ich suche speziell nach einer Möglichkeit, dies in der Datenbank zu tun und nicht durch das Schreiben eines externen Skripts.

Underdunkel
quelle
Ich denke, dass Sie dafür eine PL / Python-Funktion schreiben müssen.
Pablo
1
Hier ist ein Ausschnitt von postgis in Aktion, der helfen könnte: bostongis.com/postgis_translate.snippet
Pablo
@Pablo: Ja, höchstwahrscheinlich. Ich werde meine Frage anpassen.
underdark

Antworten:

13

Hallo

Wenn Ihre ursprüngliche Tabelle gps_p heißt, heißt Ihr Zeitstempelfeld ts und die Punkte 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;

Was es tut, ist, dass es Linien zwischen den Punkten baut und st_segmentize verwendet, um die Linie in 5 Segmente zu teilen.

Wenn zwischen den ursprünglichen Punkten nicht genau 5 Sekunden liegen, funktioniert dies nicht. Dann können Sie einfach ein ID-Feld mit einer Sequenz hinzufügen und dieses verwenden, um sich stattdessen mit id1 + 1 = id2 selbst in die Tabelle einzutragen.

HTH

/ Nicklas

Nicklas Avén
quelle
6

Hier ist ein Code-Entwurf für pl / Python. Es ist nur die Grundidee, die Punkte um eine bestimmte Distanz und einen bestimmten Azimut zu übersetzen.
Um Postgis-Funktionen in pl / python auszuführen, war die einzige Lösung, die ich gefunden habe, die Verwendung von plpy.prepare und plpy.execute (sehr langweilig).

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)
Pablo
quelle
0

Wenn ich mich nicht irre ...
Was Sie tun müssten, ist, die Verbindungslinie zu bestimmen und sie dann zu teilen.

Brad Nesom
quelle