PostGIS ST_Intersection with interpolated M value

Question

I have a polygon and a GPS Trace (incl timestamps in the M dimension) zig-zagging across the polygon.

I want to find the time intervals whenever the GPS trace entered and exited the polygon, which can be multiple times. I am doing something like:

with snippets as (
  select
    ST_Dump(ST_Intersection(t.gps_trace, u.polygon)) as snippet
  from
    table1 as t,
    table2 as u
  where
    ST_Intersects(t.gps_trace, u.polygon)
),
select
   ST_StartPoint(snippet) as enter,
   ST_EndPoint(snippet) as exit
from
  snippets
;

The problem is, that the two points enter and exit don't have their M-value retained/interpolated. Therefore, I had the idea to do something like ST_InterpolatePoint(gps_trace, ST_LocatePoint(gps_trace, enter)) to recover an interpolated M value. However, this does not always work because ST_LocatePoint looks for the spatially closest point, and if two GPS snippets enter the polygon at the same point, this will interpolate the wrong M. I also tried to use ST_SharedPaths to overlay the snippet with the original path, but that doesn't retain M either it seems (and I tend to get empty shared paths, which is weird too).

Therefore my question: is there a way to do an intersection of a linestring with a polygon that retains/interpolates M values? Or how else can I find the GPS timestamps when the line enters the polygon?

Answer 1

After trying various approaches, this is the best solution I was able to come up with. Because I saw several similar questions around the internet ("How to intersect GPS traces with polygons"), I'll post my solution here.

As mentioned in my comments, there are three main issues for a full solution:

1. Many of the PostGIS functions swallow the M-value;
2. Self-intersecting lines screw up ST_Intersection() and ST_Split(), creating spurious intersection points if the self-intersection happens inside the polygon;
3. If a GPS trace contains points with the same (x,y), for example because points are snapped to a grid, but different M values, ST_LineLocatePoint() and ST_InterpolatePoint() will return points with the wrong M value.

I was able to work around the first and second issues, the third one remains. In practice, I will simply not snap to a grid and hope that in reality it will never happen that a trace intersects the polygon boundary at the very exact point.

Here is my solution:

with
t1 as
(
SELECT 
    -- pretending this is part of something larger
    1 as id,
    -- a weird wiggly line with some self-intersections inside the polygon
    'LINESTRINGM(
        -10 2 2,0 2 2.5, 1 3 2.6, -1 3 2.7, -1 1.5 2.8, 1 2 2.9, 10 2 4,
        -8 5 6, -3 4 7, -2 3 8, -2 4 9, -6 2 10
    )'::geometry as line
    -- a triangle
    , ST_MakePolygon(ST_GeomFromText('LINESTRING(0 0,-5 5,5 5,0 0)')) as poly
),
t2 as
(
    select 
        id,
        line,
        poly,
        -- intersect line with poly boundary (not with full poly to avoid spurious
        -- intersection points!), then dump the intersection points and locate
        -- them on the original line
        ST_LineLocatePoint(
            line
          , (ST_Dump(
                ST_Intersection(line, ST_Boundary(poly))
             )).geom) as fraction
    from
        t1
),
t3 as
(
    -- t2 contains one row with the fraction for each line/poly intersection point
    -- we now pair up always two consecutive fractions, which will allow us to
    -- recreate the line segments between the intersection points.
    -- The first SELECT pairs up start/end fractions starting at 0 up 
    -- to the last fraction.
    -- This is done through the lag() window function to get the previous fraction
    -- (and a 0.0 for the first row, which doesn't have a previous).
    -- The second SELECT after the UNION adds the very last start/end pair, 
    -- from the max fraction to 1.0.
    -- In each SELECT we use those two fractions immediately 
    -- with `ST_LineSubstring()` to get the 
    -- corresponding line segment. This retains the M-value.
    select 
        id,
        line,
        poly,
        -- get the substring between two fractions
        ST_GeometryN(
            ST_LineSubstring(
                line, 
                lag(fraction, 1, 0.0) over (partition by id order by fraction), --f1
                fraction --f2
            )
            , 1
         ) as snippet
    from 
        t2
    union
    select
        id,
        line,
        poly,
        ST_GeometryN(
            ST_LineSubstring(
                line, 
                max(fraction), --f1
                1.0 --f2
            )
            , 1
        ) as snippet
    from
        t2
    group by
        id, line, poly
)
    -- now filter out the line segments that fall outside poly and we are done. 
    -- Just display the result
    select
        ST_AsText(snippet)
        , ST_Collect(ST_Buffer(snippet, 0.1), poly)
    from
        t3
    where
        -- use a buffer here, because the endpoints are exactly on top of the
        -- boundary and ST_Within and ST_Covers don't have any tolerance and
        -- may filter out wrong segments
        ST_Within(snippet, ST_Buffer(poly, 0.001))
    ;