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

[Paragon Corporation]

Simon,

This looks interesting.  Actually hadn't thought of using ST_DumpRings, but
I think that would be better than the ST_InteriorRingN

Couple of comments
1)  You should do this

SELECT (ST_DumpRings(a.geom)).*

Instead of this
SELECT (ST_DumpRings(a.geom)).geom As the_geom, path(ST_DumpRings(a.geom))
as path

(Which would mean in the upper part you would need to reference by .geom
instead of the_geom

The reason for that is internally PostgreSQL will call ST_DumpRings twice.
This was pointed out to me by a very experienced PostgreSQL developer.

His blog entry about it is here
http://www.depesz.com/index.php/2008/11/03/waiting-for-84-pl-srf-functions-i
n-selects/

2) I think ST_BuildArea might be better than ST_MakePolygon in this regard.
It will work fine with a single closed ring and if multiple, it turns the
inners to holes.

So what I was thinking in verbiage

 ST_BuildArea(ST_Collect all exterior/interior excluding all interior rings
where area < desired (that would exclude holes that are too small))  


Thanks,
Regina

-----Original Message-----
From: postgis-users-bounces at postgis.refractions.net
[mailto:postgis-users-bounces at postgis.refractions.net] On Behalf Of Simon
Greener
Sent: Thursday, January 01, 2009 1:48 AM
To: PostGIS Users Discussion
Subject: Re: [postgis-users] Help with an area and hole filter query?

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
_______________________________________________
postgis-users mailing list
postgis-users at postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users