[postgis-users] postgis-users Digest, Vol 143, Issue 15

[John Abraham]

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
jea at hbaspecto.com

On Thu, Jan 16, 2014 at 1:00 PM, <postgis-users-request at lists.osgeo.org>
wrote:Message: 3
Date: Thu, 16 Jan 2014 14:42:29 +1000
From: Paul Pfeiffer <nightdrift at gmail.com>
To: postgis-users at lists.osgeo.org
Subject: [postgis-users] Updated cleangeometry function
Message-ID:
        <CABfX3NcT=qMaw42WtZ-PUD6ZUnr7o8BMR11JB8zeFOke4Nrhmg at 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
jea at hbaspecto.com
403-232-1060
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20140917/017d09eb/attachment.html>