Splitting lines by multiple polygons, taking length of particular segment, using PostGIS

Question

I have a table of 100k+ non-overlapping POLYGON features, and a table of 50k+ simple LINESTRING features (each line has only 2 vertices). The below is a common arrangement of these features on a map, showing 1 line (red, going left to right) and 4 individual polygons (blue):

enter image description here

Imagine the above scenario, but in thousands of different locations. The lines and polygons do vary in length/size across the map (the above is just a simplified example).

I can split the line features by the polygon features by using ST_INTERSECTION(lines, polygons), so that the above example would have the line split into 4 segments.

After successfully splitting the line into 4, I want to obtain the length of only the second segment (this is the important measurement for my particular project).

I can ST_LENGTH the result of question 1, but how can I extract only the second measurement?

This is especially tricky as the ST_INTERSECTION does't output the line segments in order (left to right here). Perhaps I could somehow add a column containing an identifier for each segment, with 'a' being the first, 'b' being the second (and the one I want to measure!), then 'c', 'd' and so on... Then: SELECT ST_LENGTH FROM line_segments WHERE label = 'b'. But I need to get the segments in order somehow first...

If you want to test, here's the geometries for the example above.

Polygons (the ids are random in the real dataset):

SELECT 'ogb1' as id, ST_GeomFromText('POLYGON((2 2, 2 4, 4 4, 4 2, 2 2))') as geom
UNION
SELECT 'ogb2' as id, ST_GeomFromText('POLYGON((4 2, 4 4, 5 4, 5 2, 4 2))') as geom
UNION
SELECT 'ogb3' as id, ST_GeomFromText('POLYGON((5 2, 5 4, 5.5 4, 5.5 2, 5 2))') as geom
UNION
SELECT 'ogb4' as id, ST_GeomFromText('POLYGON((5.5 2, 5.5 4, 7 4, 7 2, 5.5 2))') as geom

and the line:

SELECT 1234 as id, ST_GeomFromText('LINESTRING(3 3, 6 3)') as geom

• dr_jts: What is the use case for this? It sounds potentially interesting, if you can share.
• Theo F: @dr_jts hello. It’s about measuring space in front of buildings to see if a particular threshold is met or not. We assume the width of this space is generally the same as the width of the building but we needed to work out the depth from pavement to building.
• Cyril Mikhalchenko: You don't happen to have ordered polygons...
• Theo F: @CyrilMikhalchenko ordered in what way? As far as I know there is no order to the polygons, but I can check to see if their id numbers are in some kind of sequential/geographical order if you like.
• Cyril Mikhalchenko: 1) what I meant was, if you knew the polygon identifier, for example 'ogb2', you would cut a line with it and then measure its length...2) so we don't have to guess, try to provide a real example of your geodata, something like that...

Answer 1

The query below will:

1 - Intersect the lines and polygons,

2 - extract the startpoints of each intersect output line geometry,

3 - Use the startpoints as input to ST_LineLocatePoint to get a number between 0-1 of where on the input lines the point is located.

4 - Order by this number and select number 2 intersection line by each original line id:

with cte as (
select lineid, row_number() over(partition by lineid order by st_linelocatepoint(linegeom, st_startpoint(intergeom))) segmentnum, st_length(intergeom) length, intergeom from 
 (select p.id polyid, 
 l.id lineid, 
 l.geom as linegeom, 
 (st_dump(st_intersection(p.geom, l.geom))).geom intergeom
 from public.poly99 p
 join public.line99 l
 on st_intersects(p.geom, l.geom)) sub
)
select * from cte
where segmentnum=2

This finishes in three seconds for 50k lines and 250k polygons:

enter image description here

( This is the code I used to generate the random 50k lines

import numpy as np
linenum = 50000 #The number of lines to create
e = iface.mapCanvas().extent() #Within current map extent
xMax, xMin, yMax, yMin = e.xMaximum(), e.xMinimum(), e.yMaximum(), e.yMinimum() #Find the areas coordinates 
vl = QgsVectorLayer("LineString?crs=EPSG:3006&index=yes", "randomLines", "memory") #Create a temp line layer. Replace 3006 with a suitable crs
provider = vl.dataProvider()
def randomline(xmin, xmax, ymin, ymax):
 """A function to return a random line within an area"""
 global point1, point2
 randx = np.random.uniform(xmin, xmax) #Create one random x coordinate
 randy = np.random.uniform(ymin, ymax) #And y
 point1 = QgsPoint(randx, randy) #Create a point from them
 randx+=np.random.uniform(-300,300) #Create a second x coordinate within +/-300 meters from the first
 randy+=np.random.uniform(-300,300) #And y
 point2 = QgsPoint(randx, randy) #And a second point
 return QgsGeometry(QgsLineString([point1,point2])) #Use them to return a line geometry
for x in range(0,linenum):
 f = QgsFeature()
 f.setGeometry(randomline(xMin, xMax, yMin, yMax))
 provider.addFeature(f)
QgsProject.instance().addMapLayer(vl)

)

• Theo F: @BERA hmm that case won’t happen as the line is always straight and the polygons fairly regular in shape. Regardless, your answer works well! Ive written it differently though, using nested sub queries instead of CTEs (which I’m not as familiar with). The test will be to apply this to huge areas of millions of polygons.
• Cyril Mikhalchenko: @BERA, If you have a lot of geodata try not to use nested multilevel queries with WITH construct, they are very costly and inefficient, as a rule on the third level or more such an SQL-construction will simply stop behaving correctly because of its increasing complexity, I just experimented with this construct in my answers and resorted to it in case of working with a small number of geobjects...
• dr_jts: It's better to use a non-nested CTE for syntactic clarity. Postgres optimizes queries across the entire CTE, so there should not be any difference in performance. You could try doing EXPLAIN ANALYZE and see.
• Theo F: @BERA answer accepted as the methodology works for me. Thank you. Also your latest image is incredible (a piece of art!). I'd love to know how you generated the random lines too...
• Bera: Nice :)! I've added the code I used to create the lines

Answer 2

ST_Intersection does keep the inherent (vertex) order of the LINESTRING, if given a compound operational context; thus, with a collection of POLYGONs to operate on, you can simply extract segments in order of the (ST_Dump).path:

SELECT
 ln.id,
 its.path[1] AS seq,
 ST_Length(its.geom) AS len,
 its.geom
FROM
 lines AS ln
 CROSS JOIN LATERAL (
 SELECT
 (ST_Dump(ST_Intersection(ln.geom, ST_Collect(pl.geom)))).*
 FROM
 polys AS pl
 WHERE
 ST_Intersects(ln.geom, pl.geom)
 ) AS its
WHERE
 its.path[1] = 2
;

• geozelot: BTW this is not meant as POC - I expect this to be the most simple and fastest solution for your task...
• dr_jts: I don't think this is correct. ST_Collect does not order the polygons along the linestring.
• Theo F: @geozelot sorry, 'POC'? Also thanks, I'll give your method a try this week and get back to you.
• geozelot: Proof Of Concept - @dr_jts pretty sure your OverlayNG keeps (linear) references while building the graph and index, and then retrieval - with preference of the first parameter? So with a collection passed in as the second, it natively returns the correct order. We had this discussion before ,) But this query actually sucks for many overlaying polygons, but outperforms other approaches for smaller counts in general.
• dr_jts: Yes, but there's no inherent spatial order to the arguments to ST_Collect.