11

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

enter image description here

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?

2
  • 1
    Are those actual inner rings of the simple polygon components of those multipolygons, or are those gaps between polygons of a multipolygon?
    – geozelot
    May 19 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
    – Cupain
    May 19 at 16:07

4 Answers 4

13

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;

Data and holes removed

3
  • 2
    It works perfectly, thanks for the very clear illustration.
    – Cupain
    May 20 at 10:04
  • 1
    Yes, it's nicer to have this bundled up as a function.
    – dr_jts
    May 22 at 5:31
  • 1
    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.
    – dr_jts
    May 22 at 5:32
7

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

13
  • 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
    – Cupain
    May 20 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.
    – geozelot
    May 20 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)
    – Cupain
    May 20 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.
    – geozelot
    May 20 at 11:18
  • 2
    @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.
    – geozelot
    May 20 at 15:28
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;
1
  • Actually, it deletes large parts of the original multipolygon. I tried dumping it first, with no better success
    – Cupain
    May 20 at 9:07
3

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

1
  • 1
    This is good, thanks!
    – dr_jts
    May 22 at 16:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.