[postgis-users] postgis-users Digest, Vol 143, Issue 15
John Abraham
jea at hbaspecto.com
Wed Sep 17 14:25:11 PDT 2014
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>
More information about the postgis-users
mailing list