[GIS] Polygon intersection with LineString in PostGis

postgispostgresqltopology

I have a table which contains GPS coordinates at certain times. I created a LINESTRING from them ordered from the oldest to newest.
I also have a POLYGON and I would like to determine the subset of my LINESTRING which lies inside my polygon. This subset should be only the points from the table – no interpolated points from the linestring. All my geometries have SRID 4326.

I used the ST_Intersection(myLinestring,myPolygon). Because the linestring goes in and out of the polygon multiple times I expect the intersection will be GEOMETRYCOLLECTION of multiple LINESTRINGs (which correspond to the individual continuous components of the path inside my polygon).

As a quick check to determine the number of components inside I ran: ST_NumGeometries(ST_Intersection(myLinestring,myPolygon)) but I got some really big number which actually exceeds even the number of coordinates in my table.
When I ran ST_AsText(ST_Intersection(myLinestring,myPolygon)) I saw that the result is rather a collection of many individual points which lie inside the polygon AND the linestring, but they ARE NOT the samples from my table.

How do I do this right?

EDIT: I need to know the connected path components lying inside the polygon, not just the isolated points inside the polygon.

I need to know the inner linestrings so I can determine the total approximate time spent inside the polygon (the time differences of the first and last point of the inner path component summed over all those components would give a precise estimate of the time spent inside the polygon).

Best Answer

If I understand your question correctly, I think @gcarillo has the right idea in his comment. Rather than building the linestrings before the intersection, determine the points inside the polygon first and then build the linestrings based on sequence groups.

Although it looks a little complex, the following query will achieve this. I've made some assumptions on the likely structure of your gps points table.

-- Build complete lines within the polygon
SELECT startTime, 
       endTime, 
       ST_MakeLine(geom_array) AS line
FROM (
    -- Aggregate the points based on sequence group
    SELECT min(captureTime) AS startTime, 
           max(CaptureTime) AS endTime, 
           array_agg(location) AS geom_array, 
           grp
    FROM (
        -- Intersect the points with a polygon and determine a group for the sequence
        SELECT *, 
               seq - ROW_NUMBER() OVER (ORDER BY captureTime) AS grp
        FROM (
            -- Sequence the points based on time captured
            SELECT captureTime, 
                   location, 
                   ROW_NUMBER() OVER (ORDER BY captureTime) AS seq 
            FROM gpsPoints
            ) orderPoints
        WHERE ST_Intersects(location,ST_MakeEnvelope(1,1,4,4)) --Replace envelope with polygon
        ) pointsInside
    GROUP BY grp
    ) pointGroup;

This will pick up points on the boundary as well as inside. Single unconnected points will be turned into a single point linestring and should be handled or removed from the query results.

Edit: Changed the sort column for the grp row_number() to get a better execution plan (one less sort).

Related Question