Finding Rings in Polygon Using Postgis

Khairu Aqsara Sudirman

Khairu Aqsara Sudirman

Dec 25, 2019 — 5 mins read
Photo by <a href="https://unsplash.com/@marjan_blan" target="_blank">Марьян Блан | @marjanblan</a> on <a href="https://unsplash.com" target="_blank">Unsplash</a>

Photo by Марьян Блан | @marjanblan on Unsplash

I have a polygon table that has many small areas and holes. Now, I would like to remove small areas and holes,I decided to have a bash at this as part of my continuing to learn PostGIS. This article will firstly look at a simple (single) polygon (I will use WKT to construct a polygon with two inner holes) which will make it easy to replicate. Then it will consider how to handle situations where a supplied polygon has no holes (inner rings) and how to handle polygons and multipolygons together before concluding with some performance figures.

 Single Polygon with Inner Rings (holes)

The function that can help us do what we want is the ST_DumpRings(geometry) function. Here we can see how the function, applied to our test geomety, will extract all the rings for processing:

SELECT ST_AsText(b.the_geom) AS final_geom, ST_Area(b.the_geom) AS area
FROM (SELECT (ST_DumpRings(a.geom)).geom AS the_geom
  FROM (SELECT ST_PolyFromText('POLYGON((0 0,20 0,20 20,0 20,0 0),(10 10,10 11,11 11,11 10,10 10),(5 5,5 7,7 7,7 5,5 5))') AS geom) a) b;

My first attempt to use ST_DumpRings to achieve the required outcome resulted in the following SQL (note how we can apply our area filter to the extracted rings):

SELECT ST_AsText(c.the_geom) AS final_geom
FROM (SELECT ST_Collect(b.the_geom) AS the_geom
FROM (SELECT (ST_DumpRings(a.geom)).geom AS the_geom
      FROM (SELECT ST_PolyFromText('POLYGON((0 0,20 0,20 20,0 20,0 0),(10 10,10 11,11 11,11 10,10 10),(5 5,5 7,7 7,7 5,5 5))') AS geom) a) b
WHERE ST_Area(b.the_geom) > 1) c;

Note that the output is a MULTIPOLYGON with each ring being a separate polygon in the collection and not a single POLYGON (as was the source data) with inner rings which is what we do want as the result.

Searching the PostGIS documentation indicated that I should be able to reconstitute the rings into a single polygon via use of the ST_Difference aggregate. Here we want to take the outer ring as the first argument to ST_Difference and subtract the set of all inner rings from it:

SELECT ST_AsEWKT(d.final_geom)
FROM 
(SELECT ST_Difference(ST_GeometryN(c.the_geom,1),
 (SELECT
 ST_Collect(ST_GeometryN(c.the_geom,
generate_series(2,ST_NumGeometries(c.the_geom)))))::geometry
) AS final_geom FROM (SELECT ST_Collect(b.the_geom) AS the_geom
	FROM (SELECT (ST_DumpRings(a.geom)).geom AS the_geom
		FROM (SELECT ST_PolyFromText('POLYGON((0 0,20 0,20 20,0 20,0 0),(10 10,10 11,11 11,11 10,10 10),(5 5,5 7,7 7,7 5,5 5))') AS geom ) a
) b WHERE ST_Area(b.the_geom) > 1) c GROUP BY c.the_geom ) d;

ERROR: set-valued FUNCTION called IN context that cannot accept a SET

This is because the ST_Collect returns geometry[] (ie a multi-geometry) and not a single geometry. (See documentation note: Do not call with a GeometryCollection as an argument) Are there other solutions? The ST_MakePolygon function looks useful: ST_MakePolygon(linestring, [linestring[]]) Creates a Polygon formed by the given shell and array of holes. You can construct a geometry array using ST_Accum.

Though it requires linestring and not polygon arguments which will result in more complex SQL (as is shown below):

SELECT ST_AsEWKT(ST_MakePolygon(c.outer_ring, d.inner_rings )) AS final_geom
FROM (/* Get outer ring of polygon */
SELECT ST_ExteriorRing(b.the_geom) AS outer_ring
  FROM (SELECT (ST_DumpRings(a.geom)).geom AS the_geom, path(ST_DumpRings(a.geom)) AS path
          FROM (SELECT ST_PolyFromText('POLYGON((0 0 1,20 0 1,20 20 1,0 20 1,0 0 1),(10 10 2,10 11 2,11 11 2,11 10 2,10 10 2),(5 5 3,5 7 3,7 7 3,7 5 3,5 5 3))') AS geom) a) b
  WHERE path[1] = 0 /* ie the outer ring */
) c,
(/* Get all inner rings > a particular area */
SELECT ST_Accum(ST_ExteriorRing(b.the_geom)) AS inner_rings
  FROM (SELECT (ST_DumpRings(a.geom)).geom AS the_geom, path(ST_DumpRings(a.geom)) AS path
          FROM (SELECT ST_PolyFromText('POLYGON((0 0 1,20 0 1,20 20 1,0 20 1,0 0 1),(10 10 2,10 11 2,11 11 2,11 10 2,10 10 2),(5 5 3,5 7 3,7 7 3,7 5 3,5 5 3))') AS geom) a) b
  WHERE path[1] > 0 /* ie not the outer ring */
    AND ST_Area(b.the_geom) > 1
) d;

Note, that we get our desired result. The splitting of the SQL into the two halves to extract the outer ring and the inner rings (separately) had to be done because the only method of reconstructing the polygon was via the ST_MakePolygon function and it needs two inputs: a linestring for the outer ring and an array of linestrings for the inner rings left after being filtered by area. Sadly, ST_Collect() on the straight output of ST_DumpRings (with no filtering by path only area) just generates a multipolygon. I tried playing around with ST_Intersection etc but these still return a multipolygon.

Now, this SQL is fine for a single, hardcoded test geometry, but I think it would be messy if we tried to use it to process a lot of polygons in a table. The best way to approach this is to encapsulate the above SQL (a complete algorithm) in a function and then use the function when processing our tabular data.

CREATE OR REPLACE FUNCTION Filter_Rings(geometry,FLOAT)
RETURNS geometry AS
$$
SELECT ST_MakePolygon(c.outer_ring, d.inner_rings) AS final_geom
FROM (/* Get outer ring of polygon */
    SELECT ST_ExteriorRing(b.the_geom) AS outer_ring
      FROM (SELECT (ST_DumpRings($1)).geom AS the_geom, path(ST_DumpRings($1)) AS path) b
      WHERE b.path[1] = 0 /* ie the outer ring */) c,
   (/* Get all inner rings > a particular area */
    SELECT ST_Accum(ST_ExteriorRing(b.the_geom)) AS inner_rings
      FROM (SELECT (ST_DumpRings($1)).geom AS the_geom, path(ST_DumpRings($1)) AS path) b
      WHERE b.path[1] > 0 /* ie not the outer ring */
        AND ST_Area(b.the_geom) > $2) d
$$ 
LANGUAGE 'sql' IMMUTABLE;

And example of the use of Filter_Rings is:

SELECT ST_AsEWKT(Filter_Rings(ST_PolyFromText('POLYGON((0 0 1,20 0 1,20 20 1,0 20 1,0 0 1),(10 10 2,10 11 2,11 11 2,11 10 2,10 10 2),(5 5 3,5 7 3,7 7 3,7 5 3,5 5 3))') ,1::FLOAT));


postgis gis
You might enjoy

Find Topology Error With PostGIS

Topology errors can be fixed quickly using the Fix Topology Error tool, This tool allows you to select a topology error and choose from a nu...