4

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:

planned track

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

correct vertical profile

Bug Description

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

intersecting linestring wrong vertical profile

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
  1. 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;
  1. 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?

3

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
8
  • I have tried your query with the simple track on the first screenshot and for the LKD10 airspace as shown in the question it returns four rows. Three have contact: false. There should be contact: true 2 times, right?
    – Andrej S
    Sep 27 at 12:20
  • @AndrejS yes. This is likely due to floating point precision issues, as I mentioned. You can use ST_DWithin instead - I'll update the query.
    – geozelot
    Sep 27 at 12:36
  • I will post a new answer here with the query I am now using. I added ST_Distance to see the distance between the point and the polygon.
    – Andrej S
    Sep 27 at 12:42
  • @AndrejS better use ST_DWithin as in my updated query.
    – geozelot
    Sep 27 at 12:42
  • 1
    @AndrejS Well, you know what? If you think about it, we can simply assume that these crossings intersect, so there's no need to compute that. Also, we don't need ST_DWithin for the start and end points. So I will add my final update.
    – geozelot
    Sep 27 at 13:51
1

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
;

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.