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 Spike-Remover function deagrades the polygon at inner and outer rings. After then the Spike Remover calls the Spike-Remover-Core. 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;
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 polygon -- 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;
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))