Top 5 Recent Articles
ARTICLES CATEGORIES
- Algorithms (22)
- All (399)
- Biography (1)
- Blog (44)
- Business Requirements (1)
- Commentary (1)
- Conversion (2)
- Customers (2)
- Data Models (1)
- Education (2)
- GeoRaptor (13)
- GPS (1)
- Image Processing (2)
- Import Export (8)
- Licensing (2)
- LiDAR (1)
- Linear Referencing (4)
- Manifold GIS (3)
- Mapping (1)
- MySQL Spatial (7)
- Networking and Routing (including Optimization) (5)
- Open Source (18)
- Oracle Spatial and Locator (194)
- Partitioning (1)
- PostGIS (36)
- Projections (1)
- Published Articles (1)
- qGIS (1)
- Recommendations (1)
- Services (1)
- Software Change Log (1)
- Source Code (37)
- Space Curves (9)
- Spatial Database Functions (109)
- Spatial DB comparison (1)
- Spatial XML Processing (11)
- SQL Server Spatial (92)
- Standards (3)
- Stored Procedure (17)
- Tessellation or Gridding (10)
- Tools (2)
- Topological Relationships (1)
- Training (2)
- Triangulation (2)
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).
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).
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.
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.
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.
Finally, what happens if we give ST_BuildArea the linework for two disjoint polygons? Here is the 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:
I hope this article was useful to someone.
Documentation
- MySQL Spatial General Functions
- Oracle LRS Objects
- Oracle Spatial Exporter (Java + pl/SQL)
- Oracle Spatial Object Functions
- Oracle Spatial Object Functions (Multi Page)
- PostGIS pl/pgSQL Functions
- SC4O Oracle Java Topology Suite (Java + pl/SQL)
- SQL Server Spatial General TSQL Functions
- SQL Server Spatial LRS TSQL Functions