SpatialDB Advisor: Building polygons from overlapping linestrings requiring intersection

Building polygons from overlapping linestrings requiring intersection

· Apr 11, 02:48 pm by Simon Greener

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.

Hey Simon,
Good article. Gave me a lot to think about.

For the case where you have 2 polygons that share an edge and you don’t want the edge dissolved as ST_BuildArea would do, I would use ST_Polygonize which will return a collection of 2 polygons and then dump to get the individual back. Though care would need to be taken in case of holes. Haven’t really thought it out too much — but would look like this.

-- 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;

— Regina Apr 12, 06:34 am #

Regina,

Thank you for pointing this out.

The article was, of course, about ST_BuildArea. But it is perfectly sensible to point out how ST_Polygonize is a complementary tool to the PostGIS practitioner in this area. I also looked at ST_MakePolygon(linestring, [linestring[]]) but its requirement to be able to identify and separate the linestring that make up the boundary and holes of polygons is very, very difficult when all one has is linestring spaghetti!

I have modified the article to include your excellent suggestion.

regards
Simon

— Simon Greener Apr 12, 12:12 pm #

<<Loading Point Data from a CSV File in PostGIS >>Converting Oracle Optimized Rectangles to PostGIS

Name		Remember
E-mail
http://
Message
	Textile Help