wiki:UsersWikiExamplesSpikeRemover

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))
Last modified 13 years ago Last modified on Apr 8, 2011, 1:00:09 AM
Note: See TracWiki for help on using the wiki.