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:
- Many of the PostGIS functions swallow the M-value;
- Self-intersecting lines screw up
ST_Split(), creating spurious intersection points if the self-intersection happens inside the polygon;
- 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_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:
-- pretending this is part of something larger
1 as id,
-- a weird wiggly line with some self-intersections inside the polygon
-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
-- 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
)).geom) as fraction
-- 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.
-- get the substring between two fractions
lag(fraction, 1, 0.0) over (partition by id order by fraction), --f1
) as snippet
) as snippet
id, line, poly
-- now filter out the line segments that fall outside poly and we are done.
-- Just display the result
, ST_Collect(ST_Buffer(snippet, 0.1), poly)
-- 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))