Summarise duration of points within polygon using PostGIS

Question

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 using SQL?

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.

Answer 1

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)