4

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?

6
  • 2
    You can adapt this solution
    – JGH
    Apr 26 at 17:24
  • Can you break GPS trace down into separate line segments (with interpolated M values) and then process each one separately? That should allow determining the correct M value for each intersection point on each segment.
    – dr_jts
    Apr 26 at 17:34
  • Thanks JGH. The solution you linked is very close to what I need, but it has two issues: 1. If the GPS track enters the polygon at the exact same location (eg because it's snapped to a grid), then ST_InterpolatePoint will find the wrong M measure. And 2. ST_Intersection() surprisingly cuts line segments when the line self-intersects inside the polygon, resulting in start/end points that are not on the polygon boundary (see stackoverflow.com/questions/62090829/…). Same behavior for ST_Split
    – Christoph
    Apr 27 at 14:16
  • Thanks dr_jts. Cutting the gps trace into line segments may solve my issue actually. If those are straight lines, then a line segment can only enter/exit the polygon at the same location once and not self-intersect (self-intersection is an issue with ST_Intersection(), see my comment above). The difficulty is then to re-assemble the correct line segments afterwards. I guess I can keep the segment# and then re-join consecutive segments if their start/end points still touch afterwards.
    – Christoph
    Apr 27 at 14:22
  • Argh, but of course ST_LineMerge will strip the M value again. This M-value business in PostGIS is a really extremely frustrating experience :(
    – Christoph
    Apr 27 at 14:30

1 Answer 1

2

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.