[postgis-users] Help with an area and hole filter query?
Obe, Regina
robe.dnd at cityofboston.gov
Thu Jan 1 19:31:48 PST 2009
What I was thinking of was this.
CREATE OR REPLACE FUNCTION upgis_filter_rings(geometry,float) RETURNS geometry AS
$$ SELECT ST_BuildArea(ST_Collect(a.geom)) as final_geom
FROM ST_DumpRings($1) AS a
WHERE a.path[1] = 0 OR
(a.path[1] > 0 AND ST_Area(a.geom) > $2)
$$
LANGUAGE 'sql' IMMUTABLE;
The main disadvantage aside from possibly speed over Simon's is that if you have 3-d polygons, I think the above will squash them to 2d where as his approach will support 3D.
On closer inspection -- when applying the above to the example polygon Simon provided the above returns a multipolygon and filter_rings returns a POLYGON and that is because if you do a validity check of the polygon, its invalid because the holes lie outside of the exterior ring. Buildarea just assumes everything outside is a polygon and everything inside is hole where as ST_MakePolygon takes your categorization of hole/vs shell as gospel.
So better compare
SELECT ST_AsText(upgis_filter_rings(ST_GeomFromText('POLYGON((10 10,10
20,20 20,20 10,10 10),(13 17, 13 14, 15 15.82, 13 17), (18 15, 18 14, 18 14, 18 15))')
,2));
SELECT ST_AsText(filter_rings(ST_GeomFromText('POLYGON((10 10,10
20,20 20,20 10,10 10),(13 17, 13 14, 15 15.82, 13 17), (18 15, 18 14, 18 14, 18 15))')
,2));
Which yield the same answer.
Hope that helps,
Regina
-----Original Message-----
From: postgis-users-bounces at postgis.refractions.net on behalf of Paragon Corporation
Sent: Thu 1/1/2009 6:36 AM
To: 'PostGIS Users Discussion'
Subject: RE: [postgis-users] Help with an area and hole filter query?
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
_______________________________________________
postgis-users mailing list
postgis-users at postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users
-----------------------------------------
The substance of this message, including any attachments, may be
confidential, legally privileged and/or exempt from disclosure
pursuant to Massachusetts law. It is intended
solely for the addressee. If you received this in error, please
contact the sender and delete the material from any computer.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20090101/7978a283/attachment.html>
More information about the postgis-users
mailing list