Summarise duration of points within polygon

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.

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

• 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