[postgis-users] Help with an area and hole filter query?

[Simon Greener]

Bob,

> I have a polygon table that has many small areas and holes.  Now, I would
> like to remove small areas and holes that are 2800 m^2.
>
> Any help or advice would be really appreciated.

I had a bash with a simple WKT Polygon with two inner holes. 

SELECT ST_AsText(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((10 10,10 20,20 20,20 10,10 10),(0 0,0 1,1 1,1 0,0 0),(5 5,5 7,7 7,7 5,5 5))') 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((10 10,10 20,20 20,20 10,10 10),(0 0,0 1,1 1,1 0,0 0),(5 5,5 7,7 7,7 5,5 5))') as geom
                       ) a
                ) b
          WHERE path[1] > 0 /* ie not the outer ring */
            AND ST_Area(b.the_geom) > 1
        ) d;

The splitting of the SQL into the two halfs 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. 

I think it might be wise to wrap this SQL up in a function if you are going to process a lot of polygons in a table.

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;

Example: select ST_AsText(Filter_Rings(ST_PolyFromText('POLYGON((10 10,10 20,20 20,20 10,10 10),(0 0,0 1,1 1,1 0,0 0),(5 5,5 7,7 7,7 5,5 5))') ,1::float));

"POLYGON((10 10,10 20,20 20,20 10,10 10),(5 5,5 7,7 7,7 5,5 5))"

Anyway, I hope it gives you some idea as to the direction you could take...

regards
SImon
-- 
SpatialDB Advice and Design, Solutions Architecture and Programming,
Oracle Database 10g Administrator Certified Associate; Oracle Database 10g SQL Certified Professional
Oracle Spatial, SQL Server, PostGIS, MySQL, ArcSDE, Manifold GIS, Radius Topology and Studio Specialist.
39 Cliff View Drive, Allens Rivulet, 7150, Tasmania, Australia.
Website: www.spatialdbadvisor.com
  Email: simon at spatialdbadvisor.com
  Voice: +613 9016 3910
Mobile: +61 418 396391
Skype: sggreener
Longitude: 147.20515 (147° 12' 18" E)
Latitude: -43.01530 (43° 00' 55" S)
NAC:W80CK 7SWP3