Hey Paul there are many problems with the clean geometry functions, so this contribution is very useful.
I'm trying it now, I notice it returns a GeometryCollection sometimes, so I had to use ST_CollectionExtract(cleangeometry(geom),3) to get my multipolygons back. The old function wasn't returning "null", but it was returning geometries where st_area() = null. Whatever that means. Yours seems to have fixed the problem in general, thank you. However comparing the ST_Area of the original and the cleaned, I see a few of my multi polygons have lost considerable area. In one case, a tiny infinitesimally narrow "finger" stuck out from the polygon, and a tiny section of the finger was retained but the bulk of the polygon was lost. In another case, it looks like the corner was twisted around itself, and so the twisted corner was retained and the bulk of the polygon was lost. I can manually make these edits, I think, and your algorithm is better than the original, but it's not perfect yet :) -- John Abraham [email protected] On Thu, Jan 16, 2014 at 1:00 PM, <[email protected]> wrote:Message: 3 Date: Thu, 16 Jan 2014 14:42:29 +1000 From: Paul Pfeiffer <[email protected]> To: [email protected] Subject: [postgis-users] Updated cleangeometry function Message-ID: <CABfX3NcT=qmaw42wtz-pud6zunr7o8bmr11jb8zefoke4nr...@mail.gmail.com> Content-Type: text/plain; charset="iso-8859-1" I have been having problems with the clean geometry function referenced at http://trac.osgeo.org/postgis/wiki/UsersWikiCleanPolygons returning null geometries. So I have fixed the script and I'm providing it below if anyone else wants to test it and let me know I could update the wiki. (Requires Postgres v9 +) -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- -- $Id: cleanGeometry.sql 2014-01-16 Paul Pfeiffer -- -- cleanGeometry - remove self- and ring-selfintersections from -- input Polygon geometries -- -- Copyright 2014 Paul Pfeiffer -- Version 2.0 -- contact: nightdrift at gmail dot com -- -- modified from cleanGeometry.sql 2008-04-24 from http://www.kappasys.ch -- -- This is free software; you can redistribute and/or modify it under -- the terms of the GNU General Public Licence. See the COPYING file. -- This software is without any warrenty and you use it at your own risk -- -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - CREATE OR REPLACE FUNCTION cleangeometry(geom "public"."geometry") RETURNS "public"."geometry" AS $BODY$ DECLARE inGeom ALIAS for $1; outGeom geometry; tmpLinestring geometry; sqlString text; BEGIN outGeom := NULL; -- Clean Polygons -- IF (ST_GeometryType(inGeom) = 'ST_Polygon' OR ST_GeometryType(inGeom) = 'ST_MultiPolygon') THEN -- Check if it needs fixing IF NOT ST_IsValid(inGeom) THEN sqlString := ' -- separate multipolygon into 1 polygon per row WITH split_multi (geom, poly) AS ( SELECT (ST_Dump($1)).geom, (ST_Dump($1)).path[1] -- polygon number ), -- break each polygon into linestrings split_line (geom, poly, line) AS ( SELECT ST_Boundary((ST_DumpRings(geom)).geom), poly, (ST_DumpRings(geom)).path[1] -- line number FROM split_multi ), -- get the linestrings that make up the exterior of each polygon line_exterior (geom, poly) AS ( SELECT geom, poly FROM split_line WHERE line = 0 ), -- get an array of all the linestrings that make up the interior of each polygon line_interior (geom, poly) AS ( SELECT array_agg(geom ORDER BY line), poly FROM split_line WHERE line > 0 GROUP BY poly ), -- use MakePolygon to rebuild the polygons poly_geom (geom, poly) AS ( SELECT CASE WHEN line_interior.geom IS NULL THEN ST_Buffer(ST_MakePolygon(line_exterior.geom), 0) ELSE ST_Buffer(ST_MakePolygon(line_exterior.geom, line_interior.geom), 0) END, line_exterior.poly FROM line_exterior LEFT JOIN line_interior USING (poly) ) '; IF (ST_GeometryType(inGeom) = 'ST_Polygon') THEN sqlString := sqlString || ' SELECT geom FROM poly_geom '; ELSE sqlString := sqlString || ' , -- if its a multipolygon combine the polygons back together multi_geom (geom) AS ( SELECT ST_Multi(ST_Collect(geom ORDER BY poly)) FROM poly_geom ) SELECT geom FROM multi_geom '; END IF; EXECUTE sqlString INTO outGeom USING inGeom; RETURN outGeom; ELSE RETURN inGeom; END IF; -- Clean Lines -- ELSIF (ST_GeometryType(inGeom) = 'ST_Linestring') THEN outGeom := ST_Union(ST_Multi(inGeom), ST_PointN(inGeom, 1)); RETURN outGeom; ELSIF (ST_GeometryType(inGeom) = 'ST_MultiLinestring') THEN outGeom := ST_Multi(ST_Union(ST_Multi(inGeom), ST_PointN(inGeom, 1))); RETURN outGeom; ELSE RAISE NOTICE 'The input type % is not supported',ST_GeometryType(inGeom); RETURN inGeom; END IF; END; $BODY$ LANGUAGE 'plpgsql' VOLATILE COST 100 ; -------------- next part -------------- An HTML attachment was scrubbed... URL: < http://lists.osgeo.org/pipermail/postgis-users/attachments/20140116/659a4495/attachment-0001.html > -- John Abraham [email protected] 403-232-1060
_______________________________________________ postgis-users mailing list [email protected] http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
