3

How can you get ST_IsSimple reason and detail similar to ST_IsValidDetail in PostGIS?

When validating geometries I can use these handy functions to not only check if geometry is valid, but also the detail when not:

SELECT reason(ST_IsValidDetail(geom)) AS reason, ST_AsText(location(ST_IsValidDetail(geom))) AS location FROM my_polylines WHERE ST_IsValid(geom) = false;

This returns reason and location of the problem. I would like to get similar information for geometries where ST_IsSimple(geom) = false. Are there any PostGIS functions to achieve that? Or any guide how to achieve it through combination of PostGIS functions?

From "Ensuring OpenGIS compliancy of geometries": A LINESTRING is simple if it does not pass through the same POINT twice (except for the endpoints, in which case it is referred to as a linear ring and additionally considered closed).

So talking about LINESTRING only, there is always only one reason - 'self-intersection' (including self-overlap). The goal is to find out how many and where exactly.

  • simplicity corresponds to validity for geometries with dimension < 2, and while officially defined for geometries with dimension = 2, PostGIS allows for others in its ST_IsValid family of functions. – geozelot Jul 8 at 7:29
  • @geozelot I have some decent number of lines with self-intersections (self-overlap) which pass through st_valid with true but do not pass st_issimple. Not really sure how your comment can help in this case. – Miro Jul 8 at 15:25
  • Hm, indeed, while it does check for validity in the sense of vertex count, the code does not imply that it would check for anything else. So, my comment is technically correct, but the implication that it would detect non-simple LineStrings cannot hold true. I was under the impression that, for convenience, it actually does extend to the definition of simple LineStrings... – geozelot Jul 8 at 15:40
3

As @dr_jts wrote, it is not as straightforward as it seems.

To find the self-intersections, one could also make use of the fact that an intersection point should be listed at least 3 times on a properly noded line (segment 1 ends at intersection point, segment 2 starts at intersection point and segment X starts or ends there too).

A properly noded line has a vertex at each intersection.

Let's note that a line whose start and end points are at the same place (a "loop") is considered simple, so we don't need to search for these cases.

So, one could

  1. properly node the input, which will add points at intersections AND change the input single line to a multi line
  2. extract each single line
  3. extract the start and end points of each single line
  4. count the number of occurrence of each extremity point
  5. return as intersection point the extremities occurring at least 3 times.
WITH data(id, geom) AS (VALUES
 ( 1, 'LINESTRING (0 0, 9 9, 0 9, 9 8, 0 7, 9 6, 0 5, 9 4, 0 3, 9 2, 0 1, 9 0)'::geometry )
,( 2, 'LINESTRING (0 0, 10 10, 10 0, 0 10)'::geometry )
),
segments (id, geom) AS (
    SELECT id, (ST_DUMP(st_node(geom))).geom 
    FROM data)
SELECT id, ST_ASTEXT(extremity_geom)
FROM (
    SELECT ID, ST_STARTPOINT(geom) as extremity_geom
    FROM segments
    UNION ALL
    SELECT ID, ST_ENDPOINT(geom) as extremity_geom
    FROM segments
    ) as extremities
GROUP BY ID, extremity_geom
HAVING COUNT(*) > 2
ORDER BY ID;

 id |     st_astext
----+--------------------
  1 | POINT(8.1 8.1)
  1 | POINT(6.3 6.3)
  1 | POINT(2.7 2.7)
  1 | POINT(3.375 3.375)
  1 | POINT(1.125 1.125)
  1 | POINT(7.875 7.875)
  1 | POINT(4.5 4.5)
  1 | POINT(0.9 0.9)
  1 | POINT(5.625 5.625)
  2 | POINT(5 5)
(10 rows)
| improve this answer | |
  • Thank you, unfortunately my polylines are combination of manual edits and some automated processing which both can cause any kind of issues, self-intersections without nodes included. – Miro Jul 10 at 10:51
  • :-) the included call to ST_NODE() will add the missing vertices automatically – JGH Jul 10 at 11:32
  • Thank you very much. I just finished some tests now. I figured out I must test both original geometry and st_node(geom) as there are problems which slips through both. I also decided to adjust the whole query to simplify things (hopefully not oversimplify :)). – Miro Jul 14 at 4:07
  • I decided to accept this answer as it helped me the most. I also added my own solution into next answer. Thanks again. – Miro Jul 16 at 0:13
4

Determining the locations which cause a LineString to be non-simple amounts to finding all the points where the line intersects itself. This is one of those spatial tasks which seems like it should be easier to do in PostGIS than it is. But currently there's no simple function call to determine the self-intersection points.

Instead, this can be done by extracting all the individual line segments from the line, and find the points where they intersect each other EXCEPT the trivial intersections which occur between adjacent segments.

Here's a SQL example to do this over multiple line records:

WITH data(id, geom) AS (VALUES
 ( 1, 'LINESTRING (0 0, 9 9, 0 9, 9 8, 0 7, 9 6, 0 5, 9 4, 0 3, 9 2, 0 1, 9 0)'::geometry )
,( 2, 'LINESTRING (0 0, 10 10, 10 0, 0 10)'::geometry )

),
segs AS (SELECT id, i, 
                ST_MakeLine( ST_PointN(geom, t.i - 1), ST_PointN(geom, t.i)) AS geom
  FROM data
  JOIN LATERAL (SELECT generate_series(2, ST_NumPoints(data.geom)) AS i ) AS t ON true 
)
SELECT a.id, ST_Intersection(a.geom, b.geom)
FROM segs a 
JOIN segs b ON ST_Intersects(a.geom, b.geom) 
  AND a.id = b.id  -- segments from same line
  AND a.i > b.i  -- process each pair of segments only once 
  AND abs(a.i - b.i) > 1;  -- don't intersect same or adjacent segments 
| improve this answer | |
  • Thank you for confirmation there is no simple way. Going to test your sql soon, just from looking at it I am not sure it will cover all the issues I have with these polylines, anyway thank you for surely savinge me a lot time to figure this code myself. – Miro Jul 10 at 10:48
1

Based on answers I have got I wrote following SQL which seems to catch all the issues I have with lines not meeting st_issimple check:

with all_lines(gid, geom) as (values 
                          ( 1, 'LINESTRING (0 0, 1 1, 0 0, 1 1)'::geometry ), 
                          ( 2, 'LINESTRING (2 2, 3 3, 3 2, 2 3)'::geometry ),
                          ( 3, 'LINESTRING (4 4, 6 6, 5 5, 0 5)'::geometry ),
                          ( 4, 'LINESTRING (7 7, 8 8)'::geometry )
                         ),
not_simple_geom as (select * from all_lines where st_issimple(geom) is distinct from true),
issues_a as (select * from (select gid, (st_dumppoints(geom)).geom geom from not_simple_geom) 
         as nodes_a group by gid, geom having count(*) > 1 order by gid),
issues_b as (select * from (select gid, (st_dumppoints(st_node(geom))).geom geom from not_simple_geom) 
         as nodes_b group by gid, geom having count(*) > 1 order by gid)
select gid, st_astext(geom) from 
         (select * from issues_a union select * from issues_b) 
         as issues order by gid;

In words first I filter lines to only the ones which do not meet geometry simplicity check. Then I test counts of extracted nodes from both original geometries, and geometries updated by function st_node (because sometimes I do not have nodes in intersections). And finally union issues as I need every location for every line only once.

| improve this answer | |

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.