2

I have table with points that are gps coordinates for multiple vehicles.

I have another table of multiple polygons that the vehicles travel through.

I’ve managed to write a python script that filters the vehicles by returning only the points inside the polygons. Now I need to figure out the duration of each vehicles inside of each polygon.

Unfortunately, each vehicle enters and exits each polygon multiple times so I can’t simply work out the difference between the max and min timestamps as that would also include time spent outside the polygons.

Is anyone able to point me in the right direction for how I might achieve this, either using SQL on its own or in combination with Python?

One way I’ve thought of is, instead of only returning the points in polygons is to add a Boolean column based on whether the row is an entry inside/outside. I’m not sure how to write this SQL query though.

I’d then sort the table on timestamp and iterate through it to look for when that value changes.

I’ve also thought of converting to lines and finding the intersects but then I’d lose the timestamp value, which is essential to the whole operation.

Any ideas gratefully received.

3

PostGIS has a limited, yet sufficient set of linear referencing functions (as well as trajectories) whose functionality is designed around the M ordinate of geometry coordinates.

In fact, your task is a lot easier to tackle when working with LineStringM, having a proper linear referencing. The M ordinate accepts any numeric value, and GPS timestamps can conveniently get converted to UNIX timestamps, if not already.

Assuming that your table gps has (change this according to your data):

  • GEOMETRY(POINT, 4326) in column geom
  • TIMESTAMP in column ts
  • INTEGER in track to denote which track a point belongs to

and your table poly

  • an id column which can get ordered
  • GEOMETRY(POLYGON, 4326) in a geom column

running

WITH
  trajectory AS (
    SELECT  track,
            ST_SetSRID(ST_MakeLine(ST_MakePointM(ST_X(geom), ST_Y(geom), EXTRACT(EPOCH FROM ts)) ORDER BY ts),4326) AS geom
    FROM    gps
    GROUP BY
            track
    )

SELECT  p.id,
        t.track,
        SUM( ST_InterpolatePoint(b.geom, ST_EndPoint(dmp.geom)) - ST_InterpolatePoint(b.geom, ST_StartPoint(dmp.geom)) )
FROM    poly AS p
JOIN    trajectory AS t
  ON    ST_Intersects(p.geom, t.geom),
        LATERAL ST_Dump(ST_Intersection(p.geom, t.geom)) AS dmp
GROUP BY
        p.id, t.track
;

should return the total time in seconds that corresponds to the parts of each track that intersect with each polygon.


This happens:

  • in the CTE, trajectories are created by aggregating points (having the UNIX timestamp assigned as M value via EXTRACT) over their corresponding track into LineStringMs

  • the polygon and (virtual) trajectory tables are joined on ST_Intersects, creating matching rows for each polygon and an intersecting trajectory

  • the LATERAL then creates the set of intersection geometries (i.e. line segments cut at polygon boundary) for each of the joined rows, expands them with ST_Dump into a list of records (cross joined with their respective base rows) and makes them available to the SELECT

  • for each of those intersecting line segments, the interpolated M value of the projected point on the initial trajectory are calculated, for both the start point and end point (which corresponds to entry and exit points) using ST_InterpolatePoint (not to be confused with ST_LineInterpolatePoint)

  • and finally their differences are summed up for each polygon and track id (GROUP BY)

| improve this answer | |
  • This looks like it could be just what I need. I've updated your variables to match mine but I've just got an error at EXTRACT: ERROR: function pg_catalog.date_part(unknown, numeric) does not exist LINE 4: ...ST_MakeLine(ST_MakePointM(ST_X(geom), ST_Y(geom), EXTRACT(EP... HINT: No function matches the given name and argument types. You might need to add explicit type casts. SQL state: 42883 Character: 115 – Owen May 23 at 12:05
  • @Owen the equivalent to my ts column needs to be of type TIMESTAMP[TZ]. Unknown usually refers to string input. How does your timestamp column look like, what type is it? – geozelot May 23 at 12:19
  • It's an integer in epoch format. I've removed the extract(EPOCH...) part as it looks like that's converting my timestamp to epoch format, which it's already in. Am I right? I'm also wondering if I need the X and Y as I already have geom column? – Owen May 23 at 12:26
  • Yeah, if it's UNIX seconds already, of course there's no need to extract it from a timestamp, tried to make it as generic as pissible with the timestamp extraction. We need to extract X and Y to create PointMs from your geometries. – geozelot May 23 at 13:01
  • Ah I see. Noticed your reply after I posted the below modification of your SQL. I'll put that back in and see how it goes. – Owen May 23 at 13:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.