Building polygons from overlapping linestrings requiring intersection

I received an email a few weeks ago asking:

I was wondering if you could post an article explaining how to create a polygon from overlapping lines, if this is possible.

I am new as in 1 week into exploring this product and am finding the documentation for ST_BuildArea and the like a little hard to understand at the moment.

This is an example of the sort of thing I would like to be able to do.

The correspondant provided me with an ascii image of 4 overlapping lines. Here is an image of the actual test data I used as it is the same as what he asked (except I have added a mid-point vertex in each linestring).

PostGIS ST_BuildArea: Starting Four Lines

First off, let’s build a table and populate it with linestrings.

 drop table test;

 create table test ( 
   gid serial );

 SELECT AddGeometryColumn('public', 'test', 'geom', 28355, 'GEOMETRY', 2);

 insert into test (geom)
 select geom
 from (select ST_SetSRID(ST_GeomFromText('LINESTRING( 1 0, 1 5, 1 11 )'),28355) as geom
       union all
       select ST_SetSRID(ST_GeomFromText('LINESTRING( 10 0, 10 5, 10 11 )'),28355) as geom
       union all
       select ST_SetSRID(ST_GeomFromText('LINESTRING( 0 1, 5 1, 11 1 )'),28355) as geom
       union all
       select ST_SetSRID(ST_GeomFromText('LINESTRING( 0 10, 5 10, 11 10 )'),28355) as geom) as g;

 select gid, ST_AsText(geom)
   from test;
gid integer st_astext text
1 LINESTRING (1 0,1 5,1 11)”
2 LINESTRING (10 0,10 5,10 11)”
3 LINESTRING (0 1,5 1,11 1)”
4 LINESTRING (0 10,5 10,11 10)”
6 POINT (1 5)”
7 POINT (1 11)”
1 POINT (1 0)”
9 POINT (10 5)”
10 POINT (10 11)”
2 POINT (10 0)”
12 POINT (5 1)”
13 POINT (11 1)”
3 POINT (0 1)”
15 POINT (5 10)”
16 POINT (11 10)”
4 POINT (0 10)”

The steps to build a polygon from the four linestrings is:

  • Intersect each line by all those lines that intersect it. Note that in SQL we have to be careful because we need to intersect each line against all the lines that possibly intersect it (see inline view named “e”).
  • We need to clip off that part of the lines that overlap an intersecting line. For this we do our intersection via ST_SymDifference as it will product a MULTILINESTRING containing all intersected segments.
  • We need to identify those segments in the SymDifference’d MULTILINESTRING linestrings that are part of our polygon. For this we must identify which elements have a start and and end point in the set of intersection points.
    • The set of intersection point (as MULTIPOINT) is created in the inline SQL named “f”.
    • The breaking up of the MULTILINESTRINGS into single linestrings representing the intersected linestrings is handled using “(ST_Dump(e.geom)).geom”.
    • The start/end point test is done using ST_StartPoint() and ST_EndPoint() against the MULTIPOINT containing the intersection points.

Here is the SQL that does all this.

 insert into test (geom)
 select ST_SetSRID(ST_BuildArea(ST_Collect(ST_GeomFromText(p.geom))),28355) as geom
   from (select DISTINCT ST_AsText(g.geom) as geom, g.iPoint 
 	   from (select (ST_Dump(e.geom)).geom,
		        f.ig as iPoint
		   from (SELECT /* This query is the set of each line with each and every line that intersects it */
			        ST_SymDifference((select a.geom from test a where a.gid = c.gid),c.geom) as geom
			   FROM (select a.gid as gid, ST_Collect(b.geom) as geom
				   from test a, 
				        test b
				  where a.gid <> b.gid
				    and ST_Intersects(a.geom,b.geom)
				  group by a.gid
				 ) as c
			 ) e, 
		        (select /* Collect the set of all intersection points in a single multipoint geometry */
		                ST_Collect(i.point) as ig
		 	   from ( select ST_SetSRID(ST_Point(ST_X(a.pnt),ST_Y(a.pnt)),28355) as point
				     from (select ST_Intersection(a.geom,b.geom) as pnt
					     from test a, 
					          test b
					    where a.gid <> b.gid
 					      and ST_Intersects(a.geom,b.geom)
				          ) as a
				  group by ST_X(a.pnt), ST_Y(a.pnt)
				  having count(*) > 1
				 ) i
			 ) f
		 ) g
   where /* We only want those linestrings that have an intersection point (see d3) at both ends */
         ST_Intersects(g.ipoint,ST_StartPoint(g.geom)) 
     and ST_Intersects(g.ipoint,ST_EndPoint(g.geom))
      ) as p;

Execution of this SQL gives us a lovely, square polygon with all the mid-point vertices in tact (I added them in to show that any internal vertices, not matter how many describing a complex line, are maintained).

PostGIS ST_BuildArea: Resultant -Orange - Polygon formed from 4 lines

Now, ST_BuildArea will work with more lines that just 4 straight lines forming a square. I will demonstrate by processing a hexagon (I won’t include the processing SQL as it is the same.)

 -- Test Hexagon
 drop table test;

 create table test ( 
   gid serial );
 SELECT AddGeometryColumn('public', 'test', 'geom', 28355, 'GEOMETRY', 2);

 insert into test (geom)
 values 
 (ST_SetSRID(ST_GeomFromText('LINESTRING(5 14, 13 10)'),28355)),
 (ST_SetSRID(ST_GeomFromText('LINESTRING(13 2, 5 -1)'),28355)),
 (ST_SetSRID(ST_GeomFromText('LINESTRING(-1 10, 7 14)'),28355)),
 (ST_SetSRID(ST_GeomFromText('LINESTRING(1 12, -2 5)'),28355)),
 (ST_SetSRID(ST_GeomFromText('LINESTRING(-2 7, 1 0)'),28355)),
 (ST_SetSRID(ST_GeomFromText('LINESTRING(-1 2, 7 -1)'),28355)),
 (ST_SetSRID(ST_GeomFromText('LINESTRING(11 12, 14 5)'),28355)),
 (ST_SetSRID(ST_GeomFromText('LINESTRING(11 0, 14 7)'),28355));

 select gid, ST_AsText(geom)
   from test;
gid integer st_astext text
1 LINESTRING (5 14,13 10)”
2 LINESTRING (13 2,5 -1)”
3 LINESTRING (-1 10,7 14)”
4 LINESTRING (1 12,-2 5)”
5 LINESTRING (-2 7,1 0)”
6 LINESTRING (-1 2,7 -1)”
7 LINESTRING (11 12,14 5)”
8 LINESTRING (11 0,14 7)”

The formed Hexagonal polygon is as follows.

PostGIS ST_BuildArea: Starting With 8 Lines

It is not possible have ST_BuildArea create two polygons where those polygons share a simple line.

Here is the data I used.

 -- Test 5 lines 2 polygons
 drop table test;

 create table test ( 
   gid serial );

 SELECT AddGeometryColumn('public', 'test', 'geom', 28355, 'GEOMETRY', 2);

 insert into test (geom)
 select geom
 from (select ST_SetSRID(ST_GeomFromText('LINESTRING( 1 0, 1 11 )'),28355) as geom
       union all
       select ST_SetSRID(ST_GeomFromText('LINESTRING( 10 0, 10 11 )'),28355) as geom
       union all
       select ST_SetSRID(ST_GeomFromText('LINESTRING( 5.5 0, 5.5 11 )'),28355) as geom
       union all
       select ST_SetSRID(ST_GeomFromText('LINESTRING( 0 1, 11 1 )'),28355) as geom
       union all
       select ST_SetSRID(ST_GeomFromText('LINESTRING( 0 10, 11 10 )'),28355) as geom) as g;

 select gid, ST_AsText(geom)
   from test;
gid integer st_astext text
1 LINESTRING (1 0,1 11)”
2 LINESTRING (10 0,10 11)”
3 LINESTRING (5.5 0,5.5 11)”
4 LINESTRING (0 1,11 1)”
5 LINESTRING (0 10,11 10)”

And you can see from the shading that there is only one polygon even though 5 linestrings went in to ST_BuildArea.

PostGIS BuildArea: 5 Lines 2 Polygons Do Not Make

However, even though this is an article about ST_BuildArea, ST_Polygonize will do what we want (thanks to Regina Obe for pointing this out in the comments below).

Here is a modified version of the SQL that finishes with the MULTIPOLYGON that ST_Polygonize produces being “exploded” into individual polygons for writing back to the table.

 -- ST_BuildArea doesn't produce two polygons so let's try ST_Polygonize
 insert into test (geom)
 SELECT geom
 FROM ST_Dump(
 (SELECT ST_SetSRID(ST_Polygonize(ST_GeomFromText(p.geom)),28355) 
   from (select DISTINCT ST_AsText(g.geom) as geom, g.iPoint 
	   from (select (ST_Dump(e.geom)).geom,
		        f.ig as iPoint
		   from (SELECT /* This query is the set of each line with each and every line that intersects it */
			        ST_SymDifference((select a.geom from test a where a.gid = c.gid),c.geom) as geom
			   FROM (select a.gid as gid, ST_Collect(b.geom) as geom
			 	   from test a, 
				        test b
				  where a.gid <> b.gid
				    and ST_Intersects(a.geom,b.geom)
				  group by a.gid
				 ) as c
			 ) e, 
		        (select /* Collect the set of all intersection points in a single multipoint geometry */
		                ST_Collect(i.point) as ig
			   from ( select ST_SetSRID(ST_Point(ST_X(a.pnt),ST_Y(a.pnt)),28355) as point
			 	     from (select ST_Intersection(a.geom,b.geom) as pnt
					     from test a, 
					          test b
					    where a.gid <> b.gid
 					      and ST_Intersects(a.geom,b.geom)
				          ) as a
				  group by ST_X(a.pnt), ST_Y(a.pnt)
				  having count(*) > 1
			  	 ) i
			 ) f
		 ) g
   where /* We only want those linestrings that have an intersection point (see d3) at both ends */
         ST_Intersects(g.ipoint,ST_StartPoint(g.geom)) 
     and ST_Intersects(g.ipoint,ST_EndPoint(g.geom))
      ) As p)) As final;

The following image shows that two polygons are created.

ST_Polygonize builds two areas as expected

Finally, what happens if we give ST_BuildArea the linework for two disjoint polygons? Here is the linework.

PostGIS ST_BuildArea: Multipolygon linework

Here is the code to create a table with the linework.

 -- Can we create multiple polygons?
 drop table test;

 create table test ( 
   gid serial );

 SELECT AddGeometryColumn('public', 'test', 'geom', 28355, 'GEOMETRY', 2);

 insert into test (geom)
 values
 ( ST_SetSRID(ST_GeomFromText('LINESTRING(-4 8, -4 14)'),28355)),
 ( ST_SetSRID(ST_GeomFromText('LINESTRING(-5 13, 2 13)'),28355)),
 ( ST_SetSRID(ST_GeomFromText('LINESTRING(-5 9, 2 9)'),28355)),
 ( ST_SetSRID(ST_GeomFromText('LINESTRING(1 14, 1 8)'),28355)),
 ( ST_SetSRID(ST_GeomFromText('LINESTRING(5 14, 5 8)'),28355)),
 ( ST_SetSRID(ST_GeomFromText('LINESTRING(4 9, 11 9)'),28355)),
 ( ST_SetSRID(ST_GeomFromText('LINESTRING(4 13, 11 13)'),28355)),
 ( ST_SetSRID(ST_GeomFromText('LINESTRING(10 14, 10 8)'),28355));

 select gid, ST_AsText(geom)
   from test;
gid integer st_astext text
1 LINESTRING
2 LINESTRING
3 LINESTRING
4 LINESTRING
5 LINESTRING
6 LINESTRING
7 LINESTRING
8 LINESTRING

Here is a modified version of the SQL that, as with ST_Polygonize, finishes with the MULTIPOLYGON that ST_BuildArea produces being “exploded” into individual polygons for writing back to the table.

 insert into test (geom)
 SELECT (ST_Dump(polys.geom)).geom
   FROM (select ST_SetSRID(ST_BuildArea(ST_Collect(p.geom)),28355) as geom
	   from (select DISTINCT g.geom, g.iPoint
		   from (select (ST_Dump(e.geom)).geom,
			         f.ig as iPoint
			   from (select /* Intersect each line witheach and every line that intersects it */
				        ST_SymDifference((select a.geom from test a where a.gid = c.gid),c.geom) as geom
				   FROM (select a.gid as gid, ST_Collect(b.geom) as geom
				    	   from test a, 
					        test b
					  where a.gid <> b.gid
					    and ST_Intersects(a.geom,b.geom)
					  group by a.gid
					 ) as c
				 ) e, 
			        (select ST_Collect(i.point) as ig
				   from ( select ST_SetSRID(ST_Point(ST_X(a.pnt),ST_Y(a.pnt)),28355) as point
					     from (select ST_Intersection(a.geom,b.geom) as pnt
						    from test a, 
						        test b
						  where a.gid <> b.gid
						    and ST_Intersects(a.geom,b.geom)
					        ) as a
					  group by ST_X(a.pnt), ST_Y(a.pnt)
					  having count(*) > 1
					 ) i
				 ) f
			 ) g
	   where ST_Intersects(g.ipoint,ST_StartPoint(g.geom)) /* We only want those linestrings that have an intersection point (see d3) at both ends */
	     and ST_Intersects(g.ipoint,ST_EndPoint(g.geom))
	      ) as p
      ) as polys;

And here are the polygons:

PostGIS ST_BuildArea: Resultant Polygons

I hope this article was useful to someone.