3

I have a PostGIS table of thousands of linestring segments wherein each line segment has only two nodes, start and end. All of the lines have a start and/or end node that matches another line segment in the table (they are all connected end-to-end). I'm trying to create a table of points representing line intersections where a line segments intersects another line segment in the table, but NOT intersections that are only intersections because they share a start or end node with another line in the table.

The query shown below selects all intersecting start and end nodes, but NOT the crossing segments (which is what I want)

CREATE TABLE intersection_points as
SELECT      
    ST_Intersection(a.geom, b.geom),
    Count(Distinct a.id)
FROM
    public."82n16nw_gsl_paths" as a,
    public."82n16nw_gsl_paths" as b
WHERE
    ST_Touches(a.geom, b.geom)
    AND a.id != b.id
GROUP BY
    ST_Intersection(a.geom, b.geom)
;

graphic of lines and intersection points: query result showing intersection points

the QGIS 'Line Intersections' tool selects the start/end node intersections AND the segment intersections that I'm after, but there is no way to only identify the segment-only points from this tool.

5
  • SELECT a.id, ST_Intersection(a.geom, b.geom) geom FROM gslpaths a JOIN gslpaths b ON ST_Crosses(a.geom, b.geom) WHERE a.id < b.id May 15, 2022 at 16:30
  • Worth noting that ST_Crosses won't find line pairs which overlap in a line.
    – dr_jts
    May 15, 2022 at 17:02
  • @dr_jts - The OP picture shows only one intersection of line segments :-)... May 15, 2022 at 17:08
  • Yes. But it turned out the actual data had overlapping lines.
    – dr_jts
    May 15, 2022 at 17:24
  • 1) @dr_jts, I did not download the data and did not analyze it; 2) I do not notice any feedback from OP; 3) by the way we should correct the name of the question "OP ask to select intersecting segments", while in the content of the question OP asks to select points, and then in the comments it turns out that points are PointZ, etc. :-)... May 15, 2022 at 19:35

1 Answer 1

4

To find lines which intersect in the interior of both lines, a query is:

CREATE TABLE intersection_points as
SELECT      
    ST_Intersection(a.geom, b.geom) AS intpt
FROM a 
JOIN b ON ST_Intersects(a.geom, b.geom)
WHERE a.id < b.id 
      AND NOT ST_Touches(a.geom, b.geom);

To find intersections where at most one of the lines can intersect at an endpoint, the query has to select line pairs which intersect but where their boundaries do not intersect (the DE-9IM pattern F):

CREATE TABLE intersection_points as
SELECT      
    ST_Intersection(a.geom, b.geom) AS intpt
FROM a 
JOIN b ON ST_Intersects(a.geom, b.geom)
WHERE a.id < b.id 
      AND ST_Relate(a.geom, b.geom, '****F****');
9
  • The first query seems promising! --it output only 6 rows, but i'm struggling to view them in qgis. I've added a primary key to the table, but when I look at the table in the DB Manager there are 4 features with POINT geometry and 2 with LINESTRING geometries. When I try to force a geometry type on the geom column I get "ERROR: Geometry has Z dimension but column does not". These should be PointZ, EPSG 3400 in my case. Input paths (as gpkg) in case you want to see those: dropbox.com/s/ohalh8h0lqdnmei/gslpaths.sql.gpkg?dl=0
    – jamierob
    May 15, 2022 at 14:38
  • My modified query: CREATE TABLE intersection_points as SELECT ST_Intersection(a.geom, b.geom) AS geom FROM gslpaths as a JOIN gslpaths as b ON ST_Intersects(a.geom, b.geom) WHERE a.id < b.id AND NOT ST_Touches(a.geom, b.geom);
    – jamierob
    May 15, 2022 at 14:41
  • 1
    It seems to be tricky to use QGIS and PostGIS data with Z ordinates. You might try using ST_Force2D(ST_Intersection(...)).
    – dr_jts
    May 15, 2022 at 16:20
  • 1
    Also, it's easier to work with the data if it has attributes. So use this: CREATE TABLE gsl_ints AS SELECT a.id AS aid, b.id AS bid, ST_Force2D( ST_Intersection(a.geom, b.geom)) AS geom FROM gslpaths AS a JOIN gslpaths AS b ON ST_Intersects(a.geom, b.geom) WHERE a.id < b.id AND NOT ST_Touches(a.geom, b.geom);. This works for me in QGIS.
    – dr_jts
    May 15, 2022 at 16:49
  • 1
    I can hardly believe it, but that last query solved my problem and it worked beautifully. I still couldn't load it into qgis, but it gave me the id's of the conflicting point and linestring geometries, which i can then easily select and delete. Thank you so much! on to the next puzzle in this project..
    – jamierob
    May 16, 2022 at 1:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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