Go to content Go to navigation and search

Home

Current PostGIS Blog Articles


Search

Browse

RSS / Atom

Email me

textpattern

Creative Commons License
All Blog Articles, Data Models and Free Source Code by Simon Greener, The SpatialDB Advisor is licensed under a Creative Commons Attribution-ShareAlike 3.0 Unported License.

Building polygons from overlapping linestrings requiring intersection

Saturday April 11 2009 at 03:48

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.

Creative Commons License

post this at del.icio.uspost this at Diggpost this at Technoratipost this at Redditpost this at Farkpost this at Yahoo! my webpost this at Windows Livepost this at Google Bookmarkspost this to Twitter

Comment [9]

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 · 11 April 2009, 19:34 · #

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 · 12 April 2009, 01:12 · #

Simon,

Great article. Once question: you mention clipping lines which intersect other lines. However if you use this input data:

insert into test (geom)select geom
from (select ST_SetSRID(ST_GeomFromText(‘LINESTRING’),28355) as geom union all select ST_SetSRID(ST_GeomFromText(‘LINESTRING’),28355) as geom union all select ST_SetSRID(ST_GeomFromText(‘LINESTRING’),28355) as geom union all select ST_SetSRID(ST_GeomFromText(‘LINESTRING’),28355) as geom union all select ST_SetSRID(ST_GeomFromText(‘LINESTRING’),28355) as geom) union all select ST_SetSRID(ST_GeomFromText(‘LINESTRING’),28355) as geom)
as g;

You will see that it does not quite work (at least as far as I can tell). Part of the problem seems to be the linestrings which show up as part of the intersection (as expected) giving a ST_X error, but even when excluding these the process does not seem to work as expected.

Any ideas why this may be?

Cheers,
James

— James Sewell · 14 January 2011, 04:13 · #

Thanks for the great article. I have a question going towards ST_Polygonize. How would you approach it if you wanted to run it on multiple polygons, and the resulting smaller polygons should bear the attribute (non-spatial) information of all the input geometries. So that the geometries that are the product of intersection have values of input geometries, those geometries that result in areas without overlap have blank values or 0s or some default , but still have the schema resultign from attributes of all inputs?
Cheers,
Martin

— Martin · 21 November 2011, 22:26 · #

Martin,

Not sure what you mean about the “smaller polygons”, but if it is an attribute reassignment issue then two methods spring to mind.

Let’s assume the data is planar enforced.

Method 1.
————-
1. For each polygon with attributes going into ST_Polygonize create a centroid point (say another geometry column in a WITH clause query) which shares the attributes with the polygon.
2. Do the ST_Polygonize using the component polygons/linestrings.
3. Merge attributes with output polygons via point in polygon with generated centroid.

Method 2.
————-
No centroid.
Use ST_Intersects/ST_Equals/ST_Overlaps etc to discover interactions of input/ouput polygons. For those with relationships then use, say, ST_Intersection with ST_Area to look for area overlaps greater than some defined parameter (eg 85%) and assign.

Just off the top of my head.

regards
Simon

Simon Greener · 22 November 2011, 00:49 · #

Thanks Simon. I am not sure how this would scale… Imagine a series of horizontal bands in one dataset (A), with attribute column X and, a set of vertical bands B with attribute Y, all in overlapping area. As a result, you would want a checkerboard of square polygons, with attribute field X,Y, all accross the area of interest. This is what I want to achieve. Of course, if you have a much less regular geometry, you end up wiht different geometries, not squares. Ultimately, I may want to use N tables, not just two.

— Martin · 22 November 2011, 01:20 · #

Martin,

Forget scaling, let’s be clearer about the problem and the approach.

I can vaguely visualise what you are suggesting and I think I can see a solution but I really do need to see some data to understand your problem. That’s the problem with spatial data…..

I have done something similar for a 3+ million cadastral dataset so solutions do exist that scale.

regards
Simon

Simon Greener · 22 November 2011, 02:07 · #

hi, i want to select all the multilinestring which intersect with a point, for this i have a line layer and a point layer. i want help in postgis query….Any help will be appreciated.

Regards
Ranjeet

— RANJEET KUMAR PASWAN · 16 April 2013, 04:56 · #

Hi Simon, great post! I’ve written a function which does similar logic described in this post to build polygons from imperfect linework

See ST_BuildAreaLinework in https://trac.osgeo.org/postgis/wiki/UsersWikiBuildPolygonsWithLines

— Mike Toews · 10 April 2015, 14:45 · #