Deleting small holes in polygons specifying the size with PostGIS

Question

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?

Answer 1

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;

Answer 2

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]POLYGONs.

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.

Answer 3

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;

Answer 4

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))!