Changes between Initial Version and Version 1 of UsersWikiExamplesSpikeRemover

Timestamp:

Dec 3, 2009, 6:30:03 AM (14 years ago)

Author:

nik1166

Comment:

—

UsersWikiExamplesSpikeRemover

@@ +1 @@
+= Spike-Remover =
+
+This function searches for spikes and removes them. Spikes arise by the intersection of geometries or by digitization errors. The geometry and the azimuth direction angle is the input. As result you get the clean geometriy back without any spikes. The process is the following: At first the SpikeRemover function deagrades the polygon at inner and outer rings. After then the Spike Remover calls the SpikeRemoverCore. This funktion cleans  the individual polygons of inner and outside spikes. As a result you get back  the clean geometry. 
+
+''Input Value: Polygon, Angle in float 0.01'' 
+
+''Return Value: Cleaned geometry'' 
+
+== Function spike remover ==
+{{{
+
+create or replace function spikeremover(geometry, angle double precision)
+  returns geometry as
+$body$ 
+select st_makepolygon(
+        (/*outer ring of polygon*/
+        select st_exteriorring(spikeremovercore($1,$2)) as outer_ring
+          from st_dumprings($1)where path[1] = 0 
+        ),  
+                array(/*all inner rings*/
+        select st_exteriorring(spikeremovercore($1, $2)) as inner_rings
+          from st_dumprings($1) where path[1] > 0) 
+) as geom
+$body$
+  language 'sql' immutable
+  cost 100;
+alter function spikeremover(geometry, double precision) owner to gis;
+}}}
+
+== Function spike remover core ==
+{{{
+
+-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+-- 
+-- $Id: spikeRemoverCore.sql 2009-10-01 08:00 Andreas Schmidt(andreas.schmidtATiz.bwl.de)  &  Nils Krüger(nils.kruegerATiz.bwl.de) $
+--
+-- spikeRemover - remove Spike from 
+-- input Polygon geometries, angle 
+-- http://www.izlbw.de/
+-- Copyright 2009 Informatikzentrum Landesverwaltung Baden-Württemberg (IZLBW) Germany
+-- Version 1.0
+--
+-- 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 spikeremovercore(geometry, angle double precision)
+  returns geometry as
+$body$declare
+ ingeom alias for $1;
+ angle  alias for $2;
+ lineusp geometry;
+ linenew geometry;
+ newgeom geometry;
+ testgeom varchar;
+ remove_point boolean;
+ newb boolean;
+ changed boolean;
+ point_id integer;
+ numpoints integer;
+
+begin
+        -- input geometry or rather set as default for the output 
+        newgeom := ingeom;
+    
+        -- check polygon
+        if (select st_geometrytype(ingeom)) = 'st_polygon' then
+            if (select st_numinteriorrings(ingeom)) = 0 then
+                
+                -- default value of the loop indicates if the geometry has been changed 
+                newb := true;   
+                --save the polygon boundary as a line 
+                lineusp := st_boundary(ingeom) as line;
+                -- number of tags
+                numpoints := st_numpoints(lineusp);
+                -- globale changevariable 
+                changed := false;
+
+                -- loop ( to remove several points)
+                while newb = true loop
+                    -- default values
+                    remove_point := false;
+                    newb := false;
+                    point_id := 0;
+                        
+
+                    -- the geometry passes pointwisely
+                    while (point_id <= numpoints) and (remove_point = false) loop
+                                -- the check of the angle at the current point of a spike including the special case, that it is the first point.
+                            if (select abs(pi() - abs(st_azimuth(st_pointn(lineusp, case when point_id= 1 then st_numpoints(lineusp) - 1 else point_id - 1 end), 
+                                st_pointn(lineusp, point_id)) - st_azimuth(st_pointn(lineusp, point_id), st_pointn(lineusp, point_id + 1))))) <= angle                then
+                                        
+
+                                    -- remove point
+                                    linenew := removepoint(lineusp, point_id - 1);
+
+                                    if linenew is not null then
+                                        raise notice '---> remove point %', point_id;
+                                        lineusp := linenew;
+                                        remove_point := true; 
+
+                                        -- if the first point is concerned, the last point must also be changed to close the line again.
+                                        if point_id = 1 then
+                                                linenew := st_setpoint(lineusp, numpoints - 2, st_pointn(lineusp, 1));
+                                                lineusp := linenew;
+                                        end if;
+                                    end if;
+                            end if;
+                            point_id = point_id + 1;
+
+                    end loop;
+
+                    -- remove point
+                    if remove_point = true then
+                        numpoints := st_numpoints(lineusp);
+                        newb := true;
+                        point_id := 0;
+                        changed := true;
+                    end if;
+
+                end loop;
+                --with the change it is tried to change back the new line geometry in a polygon. if this is not possible, the existing geometry is used
+                if changed = true then
+                    newgeom :=  st_buildarea(lineusp) as geom;
+
+                    -- errorhandling
+                    if newgeom is not null then
+                        raise notice 'creating new geometry!';
+                    else
+                        newgeom := ingeom;
+                        raise notice '-------------- area could not be created !!! --------------';
+                        testgeom:=st_astext(lineusp);
+                        raise notice 'geometry %', testgeom;
+                end if;
+            end if;
+        end if;
+    end if;
+        -- return value
+        return newgeom;
+ end;
+ $body$
+  language 'plpgsql' volatile
+  cost 100;
+alter function spikeremovercore(geometry, double precision) owner to gis;
+
+}}}
+
+== Example 1: ==
+{{{
+POLYGON((3480407.01 5483810.171,3480407.01 5483810.17,3480409.11 5483777.431,3480407.348 5483777.421,3480405.15 5483777.409,3480404.816 5483777.407,3480394.58 5483777.35,3480395.36 5483811.12,3480404.55 5483810.46,3480405.951 5483810.295,3480406.312 5483795.106,3480405.951 5483810.296,3480406.903 5483810.184,3480407.01 5483810.171))
+}}}
+'''which produces:'''
+{{{ 
+POLYGON((3480407.01 5483810.171,3480407.01 5483810.17,3480409.11 5483777.431,3480407.348 5483777.421,3480405.15 5483777.409,3480404.816 5483777.407,3480394.58 5483777.35,3480395.36 5483811.12,3480404.55 5483810.46,3480405.951 5483810.295,3480405.951 5483810.296,3480406.903 5483810.184,3480407.01 5483810.171))
+}}}
+
+----
+== Example 2: ==
+{{{
+POLYGON((3415632.49 5291021.49,3415632.488 5291021.494,3415632.49 5291021.494,3415628.93 5291028.28,3415642.95 5291001.56,3415651.18 5290985.86,3415659.27 5290984.61,3415644.71 5290947.81,3415629.17 5290921.83,3415621.28 5290929.72,3415640.21 5290959.43,3415625.38 5290971.41,3415627.79 5290983.94,3415629.49 5290992.19,3415630.14 5290995.36,3415625.65 5291022.5,3415632.49 5291021.49))
+}}}
+'''which produces:'''
+{{{ 
+POLYGON((3415632.49 5291021.49,3415632.488 5291021.494,3415632.49 5291021.494,3415642.95 5291001.56,3415651.18 5290985.86,3415659.27 5290984.61,3415644.71 5290947.81,3415629.17 5290921.83,3415621.28 5290929.72,3415640.21 5290959.43,3415625.38 5290971.41,3415627.79 5290983.94,3415629.49 5290992.19,3415630.14 5290995.36,3415625.65 5291022.5,3415632.49 5291021.49))
+}}}