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):

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

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:

( 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)

)

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
;