PostGIS - Non-simple Linestring and Polygon Intersection

Question

Application Description

I am using PostGIS as a database for a PHP application for pilots. They can plan their flights there and see all the airspace information along the track.

I have airspace data stored as Polygons and the planned track is a LineString.

The planned track can look like this, the airspace is the LKD10 rectangle, starting at ground level, going up to 1000 feet above the ground:

I have implemented a vertical profile of planned track with the airspace:

Bug Description

Everything works fine until the track intersects itself. Then the vertical profile goes mad:

The airspace should be shown along the whole chart and only once.

Logic description

This is how the logic looks like:

1. I select every airspace intersecting the linestring:

SELECT name, type, bound_lower, bound_lower_datum, bound_upper, bound_upper_datum 
FROM airspace 
WHERE ST_Intersects(location, ST_GeomFromText('LINESTRING(16.63089752197266 48.847773057644694, 16.571159362792972 48.87148983809234)', 4326)) AND bound_lower_datum NOT LIKE 'FL' AND type
NOT IN ('SECTOR', 'FIR') 
ORDER BY type ASC, name ASC

2. For each result I find the start and the end point of the intersection and locate it on the linestring:

SELECT 
 ST_LineLocatePoint(
 track::geometry, ST_StartPoint((intersection).geom)
 ) as start,
 ST_LineLocatePoint(
 track::geometry, ST_EndPoint((intersection).geom)
 ) as end
FROM
 (
SELECT track, ST_Dump(intersections) as intersection
FROM
 (
SELECT 
 track::geometry as track,
 ST_Intersection(
 track,
 (
SELECT location 
FROM airspace 
WHERE type='D' AND name='LKD10 DOLNI DUNAJOVICE')::geometry) as intersections
FROM ST_GeomFromText('LINESTRING(16.63089752197266 48.847773057644694, 16.571159362792972 48.87148983809234)', 4326) as track
 ) as a
 ) as b;

3. Then I get the intersection position in percentage and can plot it in the JS canvas correctly.

Question

I figured out that this is caused by non-simple linestring when the ST_Dump(intersections) returns more than one result for single intersection. Is there any way I can prevent this?

Answer 1

Get the intersection Points between the track and the ST_Boundary of the airspace:

• produce a table with all crossings (plus both the ST_StartPoint & ST_EndPoint)
• get their ST_LineLocatePoint fraction
• find if they ST_Intersects with the Polygon (or if they area ST_DWithin with a tiny threshold, to avoid floating point precision issues)

Exemplary query:

WITH
    track AS (
        SELECT ST_SetSRID('<track_wkt>'::GEOMETRY, 4326) AS geom
    )
    airspaces AS (
        SELECT airspace.id,
               airspace.geom
        FROM   <airspace_table> AS airspace
        JOIN   track
          ON   ST_Intersects(airspace.geom, track.geom)
    )
SELECT id,
       contact,
       frac,
       geom
FROM   (
    SELECT airspaces.id,
           ST_Intersects(startpoint, airspaces.geom) AS contact,
           0.0 AS frac,
           startpoint AS geom
    FROM   airspaces, track
           LATERAL ST_StartPoint(track.geom) AS startpoint
    UNION ALL
    SELECT airspaces.id,
           TRUE,
           ST_LineLocatePoint(track.geom, crossings.geom),
           crossings.geom
    FROM   airspaces, track,
           LATERAL ST_Dump(ST_Intersection(ST_Boundary(airspaces.geom), track.geom)) AS crossings
    UNION ALL
    SELECT airspaces.id,
           ST_Intersects(endpoint, airspaces.geom),
           1.0,
           endpoint AS geom
    FROM   airspaces, track,
           LATERAL ST_EndPoint(track.geom) AS endpoint
) q
ORDER BY
        1, 3
;

From the result set you can easily determine that

• whenever only the start and end points are returned (frac = 0.0 & frac = 1.0), the given track is fully within the respective `airspace
• whenever there are one or more crossings.geom in the set, starting with the contact state of the startpoint the inside/outside state changes with every crossing.geom

Answer 2

Based on geozelot's answer I am now using this query:

 WITH
    track AS (
        SELECT ST_SetSRID(?::GEOMETRY, 4326) AS geom
    ),
    airspaces AS (
        SELECT airspace.name as id,
               airspace.location as geom
        FROM   airspace AS airspace
        JOIN   track
          ON   ST_Intersects(airspace.location, track.geom)
        WHERE  type=? AND name=?
    )
SELECT id,
       contact,
       frac,
       geom
FROM   (
    SELECT airspaces.id,
           ST_Intersects(startpoint, airspaces.geom) AS contact,
           0.0 AS frac,
           startpoint AS geom
    FROM   airspaces, track,
           LATERAL ST_StartPoint(track.geom) AS startpoint
    UNION ALL
    SELECT airspaces.id,
           TRUE,
           ST_LineLocatePoint(track.geom, crossings.geom),
           crossings.geom
    FROM   airspaces, track,
           LATERAL ST_Dump(
                ST_Intersection(ST_Boundary(airspaces.geom::GEOMETRY), track.geom)
           ) AS crossings
    UNION ALL
    SELECT airspaces.id,
           ST_Intersects(endpoint, airspaces.geom),
           1.0,
           endpoint AS geom
    FROM   airspaces, track,
           LATERAL ST_EndPoint(track.geom) AS endpoint
) q
WHERE contact = TRUE
ORDER BY 1, 3
;