Sunday 25 November 2018

postgis - How to get smallest line segments from intersection + difference of multiple overlapping lines?


Inspired by this question, and this excellent answer.


I have multiple overlapping lines and want to extract the complete intersection of them. So from this sample data


SELECT id, geom INTO stackex FROM (VALUES
(1,ST_GeomFROMText('MULTILINESTRING((-71.1074566182813 42.3678123402605,-71.1071774374479 42.3678371063022,-71.1070346151051 42.3676745085758,-71.1069420616668 42.3675691384874))'))
,(2,ST_GeomFROMText('MULTILINESTRING((-71.1071774374479 42.3678371063022,-71.1070346151051 42.3676745085758,-71.1069420616668 42.3675691384874,-71.1068847478646 42.3675038904688,-71.1070265897396 42.3673507913022))'))

,(3,ST_GeomFROMText('MULTILINESTRING((-71.1070523458763 42.3679834790295,-71.1071774374479 42.3678371063022,-71.1070346151051 42.3676745085758))')))
as t(id, geom)

Initial Lines


Producing something like the image below, where each resulting line segment is the shortest possible intersection/difference of the inputs.


enter image description here


What I've tried:


A combination of ST_Intersection & ST_Difference


SELECT a.id,a.id||'x'||b.id  as join_id, ST_Intersection(a.geom, b.geom) AS segment
FROM stackex a

INNER JOIN stackex b ON ST_Overlaps(a.geom, b.geom)
WHERE b.id < a .id
UNION
SELECT a.id, a.id||'-'||b.id, ST_Difference(a.geom, b.geom) AS segment
FROM stackex a
INNER JOIN stackex b ON ST_Overlaps(a.geom, b.geom)
WHERE b.id <> a .id

Which almost produces the right set of lines (1-2,3-1,2x1,2-1) but misses 1x2x3 and produces larger lines such as 2-3 below and 2x1 (not pictured, but also not a great example, since 2x1 is a multiline in this case which can be split with ST_Dump). enter image description here


Would the solution be to set this up as a function with a loop until the resulting line segments are so small that there are no more overlaps? Is there a more elegant solution?



Update with more real data


There seems to be a problem with the aggregation part of the current answer when using more real-world and complex linestrings. Some of the resulting segments don't join back up with the original geometries in order to aggregate the ids. Heres the data:


SELECT id, geom INTO stackex2 FROM (VALUES
(1,ST_GeomFromText('LINESTRING(-71.059844 42.359287,-71.059851 42.359347,-71.059869 42.359462,-71.059891 42.359577,-71.059924 42.359691,-71.059992 42.359825,-71.060053 42.359945,-71.060132 42.360065,-71.06021 42.360174,-71.060324 42.360305,-71.060451 42.360426,-71.060581 42.360533,-71.060738 42.360638,-71.060912 42.360729,-71.061131 42.360824,-71.061316 42.3609,-71.061446 42.360944,-71.061609 42.360995,-71.061797 42.361053,-71.062012 42.361107,-71.062251 42.361149,-71.062469 42.361172,-71.062671 42.361182,-71.062881 42.361182,-71.062864 42.361306,-71.062873 42.361305,-71.063716 42.36129,-71.063829 42.361286,-71.064383 42.361273,-71.064554 42.361267,-71.064752 42.36126,-71.065046 42.361246,-71.065518 42.361243,-71.066361 42.361237,-71.066855 42.361233,-71.066972 42.361232,-71.067855 42.361222,-71.068729 42.3612,-71.069896 42.361168,-71.069996 42.361178,-71.070086 42.361206,-71.070309 42.361312,-71.070396 42.361348,-71.07047 42.361371,-71.070593 42.361399,-71.070682 42.361407,-71.070784 42.361406,-71.071159 42.361377,-71.071302 42.361365,-71.071599 42.361342,-71.071724 42.361336,-71.071857 42.361335,-71.071991 42.361342,-71.079054 42.361867,-71.079214 42.361879,-71.079781 42.36192,-71.081515 42.362062,-71.082379 42.362124,-71.082457 42.362127,-71.082596 42.362121,-71.082787 42.362105,-71.08289 42.362102,-71.082972 42.3621,-71.083065 42.362102,-71.083164 42.362107,-71.083245 42.362115,-71.083318 42.362126,-71.08339 42.362139,-71.083362 42.362082,-71.083343 42.362004,-71.083335 42.361956,-71.083299 42.361998,-71.083343 42.362004,-71.083362 42.362082,-71.08339 42.362139,-71.083318 42.362126,-71.083245 42.362115,-71.083164 42.362107,-71.083065 42.362102,-71.082972 42.3621,-71.08289 42.362102,-71.082787 42.362105,-71.082596 42.362121,-71.082457 42.362127,-71.082379 42.362124,-71.081515 42.362062,-71.079781 42.36192,-71.079214 42.361879,-71.079054 42.361867,-71.071991 42.361342,-71.071857 42.361335,-71.071724 42.361336,-71.071599 42.361342,-71.071412 42.361355,-71.071302 42.361365,-71.071159 42.361377,-71.071111 42.361324,-71.071079 42.36127,-71.071024 42.361124,-71.07098 42.361043,-71.070904 42.360956,-71.070738 42.360909,-71.070639 42.360892,-71.070533 42.360886,-71.070432 42.360888,-71.070312 42.360906,-71.070215 42.36093,-71.070131 42.360969,-71.07004 42.361024,-71.070037 42.360904)')),(2,ST_GeomFromText('LINESTRING(-71.096246 42.379862,-71.09602 42.379769,-71.095897 42.379759,-71.09579 42.379721,-71.095696 42.379711,-71.095567 42.379631,-71.094337 42.379181,-71.094215 42.379133,-71.094105 42.379098,-71.093998 42.379063,-71.093871 42.379016,-71.093568 42.378902,-71.093098 42.378718,-71.092866 42.378635,-71.092589 42.378536,-71.092457 42.37849,-71.09222 42.378412,-71.091857 42.378292,-71.091043 42.377988,-71.090633 42.377845,-71.090162 42.377677,-71.089803 42.377519,-71.089267 42.377284,-71.089231 42.37715,-71.089155 42.37687,-71.089119 42.376737,-71.089046 42.376428,-71.088948 42.376189,-71.08878 42.37582,-71.088424 42.374964,-71.08837 42.37484,-71.088206 42.374496,-71.08806 42.374156,-71.088005 42.374087,-71.087935 42.374024,-71.087647 42.373768,-71.087408 42.373597,-71.086987 42.373276,-71.086892 42.373205,-71.086461 42.373053,-71.085213 42.372581,-71.08462 42.37236,-71.0841 42.372297,-71.083721 42.372251,-71.082967 42.372147,-71.0822 42.372045,-71.081339 42.371928,-71.080508 42.371821,-71.07968 42.371713,-71.079483 42.371687,-71.078441 42.371551,-71.078348 42.371544,-71.078273 42.371549,-71.078209 42.371559,-71.078109 42.37158,-71.078051 42.371593,-71.077812 42.371524,-71.077679 42.371474,-71.077649 42.371428,-71.077609 42.371344,-71.077547 42.371298,-71.077039 42.371037,-71.07682 42.370929,-71.076696 42.370865,-71.076636 42.370833,-71.076704 42.370689,-71.076743 42.370606,-71.076923 42.370636,-71.077066 42.370027,-71.076948 42.370011,-71.076612 42.369967,-71.076434 42.369943,-71.076071 42.369895,-71.075759 42.369854,-71.075704 42.370107,-71.075604 42.370097,-71.075499 42.370076,-71.075384 42.370042,-71.075275 42.37,-71.075164 42.370125,-71.075048 42.370066,-71.075106 42.370002,-71.075339 42.369745,-71.075434 42.369644,-71.0755 42.369575,-71.074016 42.368801,-71.073904 42.368744,-71.074023 42.368623,-71.074104 42.368683,-71.07458 42.368927,-71.074761 42.369029,-71.075029 42.369173,-71.075452 42.369373,-71.075561 42.369422,-71.075656 42.369372,-71.075723 42.369405,-71.075824 42.36943,-71.07588 42.369436,-71.075981 42.369434,-71.076043 42.369424,-71.076141 42.369393,-71.076255 42.369322,-71.076335 42.369221,-71.076362 42.369108,-71.076334 42.368994,-71.076284 42.368923,-71.076222 42.368869,-71.07614 42.368823,-71.076031 42.36879,-71.075927 42.368779,-71.075777 42.368795,-71.075693 42.368823,-71.075643 42.368848,-71.07558 42.368892,-71.075522 42.368955,-71.075472 42.368997,-71.0754 42.369006,-71.074722 42.368664,-71.074738 42.368644,-71.074684 42.368616,-71.073908 42.36821,-71.073819 42.368107,-71.073285 42.36783,-71.073279 42.367777,-71.073694 42.367346,-71.073752 42.367329,-71.073824 42.367329,-71.073819 42.367382,-71.073818 42.36744,-71.073823 42.367496,-71.073835 42.367552,-71.073854 42.367604,-71.073888 42.367656,-71.073933 42.36771,-71.073998 42.367776,-71.074039 42.367828,-71.074061 42.367872,-71.074095 42.367956,-71.074108 42.368044,-71.074072 42.36812,-71.074016 42.368175,-71.073631 42.368547,-71.073473 42.368663,-71.073246 42.368884,-71.07313 42.368912,-71.073012 42.368925,-71.072884 42.368907,-71.072799 42.368884,-71.072423 42.36871,-71.071684 42.368346,-71.070729 42.367864,-71.06967 42.367319,-71.069505 42.367236,-71.069279 42.367122,-71.069057 42.367012,-71.06842 42.366678,-71.068738 42.366384,-71.068813 42.366349,-71.069353 42.365872,-71.069455 42.365784,-71.069493 42.365744,-71.069532 42.36569,-71.069597 42.365572,-71.069771 42.36521,-71.070027 42.36468,-71.070385 42.363923,-71.070428 42.363848,-71.070477 42.363778,-71.070549 42.3637,-71.070625 42.36363,-71.070702 42.363574,-71.071661 42.362957,-71.071911 42.362792,-71.072016 42.362712,-71.072111 42.362621,-71.072326 42.362366,-71.072459 42.362204,-71.072504 42.362154,-71.072619 42.362009,-71.072681 42.361914,-71.072742 42.361785,-71.072789 42.361662,-71.072824 42.361553,-71.072896 42.361145,-71.0729 42.361067,-71.072902 42.360992,-71.072896 42.36094,-71.072878 42.36089,-71.072836 42.360803,-71.072739 42.360652,-71.072774 42.360644,-71.072867 42.36062,-71.072989 42.360564,-71.073046 42.360629,-71.072927 42.360679,-71.07284 42.360716,-71.072753 42.360752,-71.072675 42.360778,-71.072594 42.360801,-71.072501 42.360825,-71.072409 42.360842,-71.072305 42.360861,-71.072174 42.360879,-71.07208 42.360885,-71.071981 42.360887,-71.071824 42.360885,-71.071608 42.360867,-71.071496 42.360881,-71.071728 42.360926,-71.071778 42.360826,-71.071563 42.360887,-71.071445 42.360903,-71.07133 42.360907,-71.071134 42.360883,-71.071145 42.36082,-71.071024 42.360804,-71.070882 42.360786,-71.070687 42.360781,-71.07052 42.360802,-71.070407 42.360848,-71.070312 42.360906,-71.070215 42.36093,-71.070131 42.360969,-71.07004 42.361024,-71.069965 42.361058,-71.069888 42.361072,-71.069808 42.361075,-71.06946 42.361076,-71.068745 42.36109,-71.068257 42.3611,-71.067856 42.361109,-71.066958 42.361118,-71.066383 42.361128,-71.066002 42.361129,-71.065515 42.361134,-71.065156 42.361136,-71.064557 42.361147,-71.064152 42.361156,-71.063816 42.361166,-71.063723 42.361168,-71.06337 42.361176,-71.062881 42.361182,-71.062671 42.361182,-71.062469 42.361172,-71.062251 42.361149,-71.062012 42.361107,-71.061797 42.361053,-71.061609 42.360995,-71.061446 42.360944,-71.061316 42.3609,-71.061131 42.360824,-71.060912 42.360729,-71.060738 42.360638,-71.060581 42.360533,-71.060451 42.360426,-71.060324 42.360305,-71.06021 42.360174,-71.060132 42.360065,-71.060053 42.359945,-71.059992 42.359825,-71.059924 42.359691,-71.059891 42.359577,-71.059869 42.359462,-71.059851 42.359347,-71.059847 42.359226,-71.059849 42.359137,-71.059859 42.359026,-71.059874 42.358941,-71.059892 42.358831,-71.059952 42.358715,-71.060021 42.358589)')),(3,ST_GeomFromText('LINESTRING(-71.069246 42.361184,-71.069896 42.361168,-71.069996 42.361178,-71.070086 42.361206,-71.070309 42.361312,-71.070396 42.361348,-71.07047 42.361371,-71.070593 42.361399,-71.070682 42.361407,-71.070784 42.361406,-71.071159 42.361377,-71.071302 42.361365,-71.071599 42.361342,-71.071724 42.361336,-71.071857 42.361335,-71.071991 42.361342,-71.079054 42.361867,-71.079214 42.361879,-71.079781 42.36192,-71.081515 42.362062,-71.082379 42.362124,-71.082457 42.362127,-71.082596 42.362121,-71.082787 42.362105,-71.08289 42.362102,-71.082972 42.3621,-71.083065 42.362102,-71.083164 42.362107,-71.083245 42.362115,-71.083318 42.362126,-71.08339 42.362139,-71.083362 42.362082,-71.083343 42.362004,-71.083299 42.361998,-71.083189 42.361983,-71.082914 42.361964,-71.082808 42.361955,-71.080547 42.361784,-71.080422 42.361774,-71.080316 42.361644,-71.080395 42.361475,-71.081293 42.361237,-71.0815 42.361175,-71.081621 42.361137,-71.081657 42.361198,-71.081939 42.361646,-71.081901 42.361714,-71.081979 42.361766,-71.081963 42.361926,-71.082486 42.361966,-71.082655 42.361998,-71.082723 42.362034,-71.082787 42.362105,-71.082596 42.362121,-71.082457 42.362127,-71.082379 42.362124,-71.081515 42.362062,-71.079781 42.36192,-71.079214 42.361879,-71.079054 42.361867,-71.074732 42.361547,-71.071991 42.361342,-71.071857 42.361335,-71.071724 42.361336,-71.071599 42.361342,-71.071302 42.361365,-71.071159 42.361377,-71.071111 42.361324,-71.071079 42.36127,-71.071024 42.361124,-71.07098 42.361043,-71.070904 42.360956,-71.070738 42.360909,-71.070639 42.360892,-71.070533 42.360886,-71.070432 42.360888,-71.070312 42.360906,-71.070215 42.36093,-71.070131 42.360969,-71.07004 42.361024,-71.070037 42.360984)'))
as a(id,geom)

Using ST_Overlaps as an additional join for segments that don't join using ST_Contains works on some, but not all of the segments that don't join using ST_Contains



Answer



You can get the smallest intersections of the linework by splitting the unioned linework geometry (blue) by the unioned boundary points from the linework (red).


result



In SQL it looks like this, with result in WKT:


SELECT ST_Split(ST_Union(geom), ST_Union(ST_Boundary(geom)))
FROM stackex;

GEOMETRYCOLLECTION(
LINESTRING(-71.1071774374479 42.3678371063022,-71.1070346151051 42.3676745085758),
LINESTRING(-71.1070346151051 42.3676745085758,-71.1069420616668 42.3675691384874),
LINESTRING(-71.1069420616668 42.3675691384874,-71.1068847478646 42.3675038904688,-71.1070265897396 42.3673507913022),
LINESTRING(-71.1074566182813 42.3678123402605,-71.1071774374479 42.3678371063022),
LINESTRING(-71.1070523458763 42.3679834790295,-71.1071774374479 42.3678371063022))




And for bonus points, if you need to label the join_id from the original IDs (as in the figure above), here is how it can be aggregated:


WITH smallest_segements AS (
SELECT (ST_Dump(ST_Split(ST_Union(geom), ST_Union(ST_Boundary(geom))))).*
FROM stackex
)
SELECT row_number() over() AS rn, join_id, s.geom AS segment
FROM smallest_segements s, LATERAL (
SELECT string_agg(id::text, '+') AS join_id

FROM stackex a
WHERE ST_Contains(a.geom, s.geom)
) l;

rn | join_id | segment
----+---------+----------------------------------------------------------------------------------------------------------------------
1 | 1+2+3 | LINESTRING(-71.1071774374479 42.3678371063022,-71.1070346151051 42.3676745085758)
2 | 1+2 | LINESTRING(-71.1070346151051 42.3676745085758,-71.1069420616668 42.3675691384874)
3 | 2 | LINESTRING(-71.1069420616668 42.3675691384874,-71.1068847478646 42.3675038904688,-71.1070265897396 42.3673507913022)
4 | 1 | LINESTRING(-71.1074566182813 42.3678123402605,-71.1071774374479 42.3678371063022)

5 | 3 | LINESTRING(-71.1070523458763 42.3679834790295,-71.1071774374479 42.3678371063022)
(5 rows)

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...