The query should use ST_UNION to make sure that the linestrings are clipped at all the points they intersect and a lateral join so it would use the proper indexes.
SELECT row_number() over() AS gid,
ST_CollectionExtract(ST_Multi(ST_Difference(a.geom, b.geom)), 2)::geometry(MultiLineString, 27700) as geom
FROM lines a
,LATERAL (
SELECT ST_UNION(polygons.geom)
FROM polygons
WHERE ST_Intersects(a.geom,polygons.geom)
) AS b
Here is what we can do,
- Create a
GeometryCollection
with one polygon using the ST_Polygonize
aggregate. This is our test data.
- Access the first element of the
GeometryCollection
with ST_GeometryN(geom,1)
. If your input is already a polygon, you'll pick up from here.
- Extract the boundary with
ST_Boundary
. The boundary is a LINESTRING
.
- Decompose that
LINESTRING
into points using ST_PointN
, and using those points as arguments to ST_MakeLine
It looks like this,
SELECT ST_AsText(ST_Boundary(ST_GeometryN(ST_Polygonize(l),1)))
FROM ( VALUES
( ST_MakeLine(ST_MakePoint(5,0),ST_MakePoint(0,0)) ),
( ST_MakeLine(ST_MakePoint(0,0),ST_MakePoint(0,5)) ),
( ST_MakeLine(ST_MakePoint(0,5),ST_MakePoint(5,5)) ),
( ST_MakeLine(ST_MakePoint(5,5),ST_MakePoint(5,0)) )
) AS t(l);
st_astext
---------------------------------
LINESTRING(5 0,0 0,0 5,5 5,5 0)
(1 row)
Now we have to decompose the boundary into 2-point line strings. First we move the aggregate into an inner-select
SELECT ST_AsText(b)
FROM (
SELECT ST_Boundary(ST_GeometryN(ST_Polygonize(l),1)) AS b
FROM ( VALUES
( ST_MakeLine(ST_MakePoint(5,0),ST_MakePoint(0,0)) ),
( ST_MakeLine(ST_MakePoint(0,0),ST_MakePoint(0,5)) ),
( ST_MakeLine(ST_MakePoint(0,5),ST_MakePoint(5,5)) ),
( ST_MakeLine(ST_MakePoint(5,5),ST_MakePoint(5,0)) )
) AS t(l)
) AS i(b);
Then we recompose simple lines,
SELECT ST_AsText(ST_MakeLine(
ST_PointN(b.geom,gs.pt),
ST_PointN(b.geom,gs.pt+1)
))
FROM (
SELECT ST_Boundary(ST_GeometryN(ST_Polygonize(l),1)) AS b
FROM ( VALUES
( ST_MakeLine(ST_MakePoint(5,0),ST_MakePoint(0,0)) ),
( ST_MakeLine(ST_MakePoint(0,0),ST_MakePoint(0,5)) ),
( ST_MakeLine(ST_MakePoint(0,5),ST_MakePoint(5,5)) ),
( ST_MakeLine(ST_MakePoint(5,5),ST_MakePoint(5,0)) )
) AS t(l)
) AS b(geom)
CROSS JOIN LATERAL generate_series(1,ST_NPoints(b.geom)-1) gs(pt);
That returns..
st_astext
---------------------
LINESTRING(5 0,0 0)
LINESTRING(0 0,0 5)
LINESTRING(0 5,5 5)
LINESTRING(5 5,5 0)
(4 rows)
Best Answer
This might be a good spot to use an SQL-language function. Here's a quick one that should work for this situation:
This will retain the polygonal components of an intersection, but throw away everything else. It always returns a MultiPolygon, even if you have one or no components.