# Deleting small holes in polygons specifying the size with PostGIS

I have a big MultiPolygon layer in which I want to delete holes smaller than 100m².

QGIS do it well and fast with the algorithm `"native:deleteholes"`, but I need to use it as part of a SQL workflow with PostGIS. I found this:

``````SELECT gid, ST_Collect(ST_MakePolygon(geom)) As geom
FROM (
SELECT gid, ST_ExteriorRing((ST_Dump(geom)).geom) As geom
FROM my_spatial_table
) s
GROUP BY gid;
``````

But it does not allow to specify the max dimension of the holes to delete.

Is there a way?

• Are those actual inner rings of the simple polygon components of those multipolygons, or are those gaps between polygons of a multipolygon? Commented May 19, 2022 at 15:49
• In this case, it's one huge and complex multipolygon (created by a ST_Union of multiples multipolygons). In this unique multipolygon, there are some gaps (but I want to keep it), inner and outer spikes (I want to keep it too), and holes. I want to delete holes smaller than 100m² to simplify the geometry without loosing too much spatial information, in order to perform a ST_Difference with another layer after that Commented May 19, 2022 at 16:07

This is one of those processes which are simple in an iterative language, but trickier in set-oriented SQL. But the `ST_DumpX` set-returning functions in PostGIS and the array handling in Postgres make it easier (but still fiddly to get working). The following solution works on MultiPolygons and also Polygons if they happen to be in the dataset.

``````WITH data(id, geom) AS (VALUES
(1, 'MULTIPOLYGON (((100 100, 100 0, 0 0, 0 100, 100 100), (10 10, 10 70, 60 10, 10 10), (30 90, 90 90, 90 30, 30 90), (20 80, 10 80, 10 90, 20 80), (90 10, 80 10, 80 20, 90 10)), ((0 170, 100 170, 100 120, 0 120, 0 170), (10 130, 10 140, 20 130, 10 130)))'::geometry)
,(2, 'MULTIPOLYGON (((200 100, 300 100, 300 0, 200 0, 200 100), (210 10, 210 70, 260 10, 210 10), (280 80, 280 90, 290 80, 280 80)), ((200 160, 260 160, 260 120, 200 120, 200 160)))'::geometry)
,(3, 'POLYGON ((110 90, 190 90, 190 10, 110 10, 110 90), (120 20, 120 80, 180 20, 120 20), (170 70, 170 80, 180 70, 170 70))'::geometry)
)
SELECT id, ST_Collect(
ARRAY( SELECT ST_MakePolygon(
ST_ExteriorRing(geom),
ARRAY( SELECT ST_ExteriorRing( rings.geom )
FROM ST_DumpRings(geom) AS rings
WHERE rings.path[1] > 0 AND ST_Area( rings.geom ) >= 100
)
)
FROM ST_Dump(geom) AS poly )
) AS geom
FROM data;
``````

• It works perfectly, thanks for the very clear illustration. Commented May 20, 2022 at 10:04
• Yes, it's nicer to have this bundled up as a function. Commented May 22, 2022 at 5:31
• It's better to post that as a separate answer, though, so that authorship is clear, and so it can be commented and voted on separately. Commented May 22, 2022 at 5:32
• How would this work with ogr2ogr and the SQLite dialect? I already tried this query and the function from the other answer but all I get is `ERROR 1: In ExecuteSQL(): sqlite3_prepare_v2(...) near "SELECT": syntax error` I need this to remove small holes from a shapefile. Commented Mar 22 at 14:07

Just to add a little less convoluted solution in SQL, shifting construction logic to `ST_BuildArea` - the `LATERAL` join predicate is an elegant method to expand and iterate over sets of rows returned from functions or statements:

``````SELECT
gid,
ST_BuildArea(ST_Collect(r_dmp.geom) FILTER( WHERE r_dmp.path[1] = 0 OR ST_Area(r_dmp.geom) > 100 )) AS geom
FROM
<polygons>,
LATERAL ST_Dump(geom) AS g_dmp,
LATERAL ST_DumpRings(g_dmp.geom) AS r_dmp
GROUP BY
gid
;
``````

Fully generic and working on all collections of (previously valid) `[MULTI]POLYGON`s.

While `ST_BuildArea` mitigates the need for fiddly SQL here, it has to apply various sanity and overlap checks in order to construct valid areal geometries from arbitrary lineworks and is thus comparably slow.

Note:

As with all geometric property and predicate functions in PostGIS, `ST_Area` operates on the unit of measurement of the underlying CRS. To be able to filter by `m²` you will either have to use a suitable projection, or cast your `GEOMETRY` types to `GEOGRAPHY` on-the-fly, i.e. `WHERE ST_Area(geom::GEOGRAPHY) > 100`.

• Indeed, very elegant solution, but it requires to dump first the multipolygon, and to union the result to get a unique object as i need Commented May 20, 2022 at 10:03
• @Cupain then your input geometry is invalid. This query dumps any collection into it's atomic parts and re-collects them based on their `gid` - and with valid input areas you will always end up with a collection of one or more (former shell) polygons containing any number of (former ring) polygons (including islands within inner rings). `ST_Union` would only ever make sense if any of these overlap (apart from the inner ring containment) - which would render an areal geometry invalid as per OGC definition. Commented May 20, 2022 at 10:32
• Indeed, you're right, no need to dump. As I have only 1 feature (so 1 gid), I had tried without the "GID" in your query. With the GID, it leads to the same result as @dr_jts metod, but not as faster (4s v/s 60s) Commented May 20, 2022 at 11:05
• @Cupain ah, just move the filter to a regular `WHERE` clause. I made that change to the query above, should have similar performance now. Commented May 20, 2022 at 11:18
• @Cupain yeah, I just had a second to test run on obscenely large MultiPolygons, and its plenty slower - and that is no wonder, `ST_BuildArea` is a very costly function. For posterity I added a high performant query, structurally somewhere in between dr_jts fiddly answer and this more elegant construct. Commented May 20, 2022 at 15:28

Based on dr_jts response, I propose the following version of the design and name of the new spatial SQL function:

Create a spatial custom SQL function:

``````CREATE OR REPLACE FUNCTION ST_RemovingHolesInPolygonsByArea(
geom GEOMETRY,
area real)
RETURNS GEOMETRY AS
\$BODY\$
WITH
tbla AS (SELECT ST_Dump(geom))
SELECT ST_Collect(ARRAY(SELECT ST_MakePolygon(ST_ExteriorRing(geom),
ARRAY(SELECT ST_ExteriorRing(rings.geom) FROM ST_DumpRings(geom) AS rings
WHERE rings.path[1]>0 AND ST_Area(rings.geom)>=area))
FROM ST_Dump(geom))) AS geom FROM tbla;
\$BODY\$
LANGUAGE SQL;
``````

RUN

``````SELECT ST_RemovingHolesInPolygonsByArea(geom, 100) geom FROM <geodata>
``````

Note: The author and developer of the function itself is Martin Davis (dr jts (https://gis.stackexchange.com/users/14766/dr-jts))!

• This is good, thanks! Commented May 22, 2022 at 16:42

Based on your example, just extract the interior ring of your geom, calculate it's area and use it as condition in where clause as:

``````SELECT gid, ST_Collect(ST_MakePolygon(geom)) As geom
FROM (
SELECT gid, ST_ExteriorRing((ST_Dump(geom)).geom) As geom,
st_area(st_makepolygon((st_interiorringn((ST_Dump(geom)).geom,1)))) as interior_ring_area
FROM my_patial_table
) s
where interior_ring_area>=100
GROUP BY gid;
``````
• Actually, it deletes large parts of the original multipolygon. I tried dumping it first, with no better success Commented May 20, 2022 at 9:07