Finding and Editing Polygon Overlaps
The function pgoverlap checks a polygon table for overlapping geometries, separating between sliver and non-sliver overlappings. It returns two tables - one table with non-overlapping geometries and adjusted sliver-overlapping polygons and one table with non-sliver overlappings which have to be handled manually afterwards.
This function is a test version. Any comments or ideas for improvement are welcome!
Preconditions for successful use of the function:
- gid column of data type integer or serial, unique, with btree index!!
- the_geom column of data type geometry and geometry type POLYGON, with gist index!!
- column of data type integer (name is variable) for subsequent joining of attributes with btree index!!
- geometries should be valid (otherwise number of exceptions rises)
Caution: Do not call the functions twice at a time!! - Intermediate tables have fixed names and are locked during function processing (otherwise their content would be overwritten).
Example:
CREATE TABLE pgoverlap_bsp_input (gid serial PRIMARY KEY, a_id integer, the_geom geometry);
INSERT INTO pgoverlap_bsp_input (a_id, the_geom) VALUES
(1,ST_GeomFromText('POLYGON((3565961 5361731,3570432 5359825,3574390 5362977,3577469 5360631,3582526 5360192,3591616 5351909,3592715 5343113,3591835 5337395,3586265 5339448,3579008 5336882,3572045 5340254,3567500 5338788,3564421 5342013,3562296 5353301,3566327 5356527,3565961 5361731))',31467)),
(2,ST_GeomFromText('POLYGON((3587846 5359101,3587353 5362736,3592036 5366371,3598813 5365755,3601708 5361997,3592528 5358855,3587846 5359101))',31467)),
(3,ST_GeomFromText('POLYGON((3600429 5353669,3606823 5355001,3612240 5361751,3622453 5355801,3637461 5355179,3635951 5340171,3617036 5337240,3605225 5340260,3600429 5353669))',31467)),
(4,ST_GeomFromText('POLYGON((3580359 5328360,3588529 5330935,3597943 5323120,3606024 5322587,3612951 5315217,3608155 5306691,3603804 5305093,3596166 5306159,3591105 5310244,3584799 5310688,3576540 5317437,3580359 5328360))',31467)),
(5,ST_GeomFromText('POLYGON((3568123 5343184,3577776 5350721,3586922 5344624,3584890 5336664,3593104 5330312,3586499 5322014,3572187 5325232,3568123 5343184))',31467)),
(6,ST_GeomFromText('POLYGON((3594628 5349620,3601741 5358935,3617153 5357580,3621387 5348519,3617153 5339712,3610040 5334293,3616137 5322945,3607923 5311598,3600302 5314392,3596068 5319727,3600725 5329974,3596406 5339882,3594628 5349620))',31467)),
(7,ST_GeomFromText('POLYGON((3622319 5325401,3625029 5327433,3630194 5329466,3635191 5326502,3638747 5317610,3631719 5311598,3622150 5316594,3622319 5325401))',31467));
SELECT _pgoverlap('public','pgoverlap_bsp_input','a_id','pgoverlap_bsp_output','pgoverlap_bsp_todo');
Function:
CREATE OR REPLACE FUNCTION _pgoverlap(schema_a varchar, table_a varchar, a_id varchar, table_result1 varchar, table_result2 varchar)
RETURNS void AS
$BODY$
/*
$Id: pgoverlap_wiki.sql 2011-03-09 15:19Z Birgit Laggner $
pgoverlap - checks table for overlapping geometries separating between sliver and non-sliver overlappings, returns table with non-overlapping geometries and adjusted sliver-overlapping polygons as well as table with non-sliver overlappings (for later treatment)
schema_a: database schema where table is located
table_a: table name
a_id: id column name in table
table_result1: name of 1st result table (non-overlapping polygons)
table_result2: name of 2nd result table (overlapping polygons)
Preconditions for successful use of the function:
- gid column of data type integer or serial, unique, with btree index!!
- the_geom column of data type geometry and geometry type POLYGON, with gist index!!
- column of data type integer (name is variable) for subsequent joining of attributes, with btree index!!
- geometries should be valid (otherwise number of exceptions rises)
Caution: Do not call the function twice at a time!! - Intermediate tables have fixed names and are locked during function processing (otherwise their content would be overwritten).
Copyright (C) 2011 Johann Heinrich von Thünen-Institute (vTI) - Federal Research Institute for Rural Areas, Forestry and Fisheries, Institute of Rural Studies, Braunschweig, Germany (http://www.vti.bund.de)
Version 0.1
contact: birgit dot laggner at vti dot bund dot de
This is free software; you can redistribute and/or modify it under
the terms of the GNU General Public Licence. This software is without any warrenty and you use it at your own risk.
*/
DECLARE
i integer;
sql_1 text;
BEGIN
--find all overlappings with geometrytype Polygon:
EXECUTE
'DROP TABLE IF EXISTS '||schema_a||'.pgoverlap_tmp1;';
EXECUTE
'CREATE TABLE '||schema_a||'.pgoverlap_tmp1 (
gid serial,
a_'||a_id||' integer,
b_'||a_id||' integer,
a_geom geometry,
b_geom geometry,
i_geom geometry,
area numeric,
perimeter2 numeric,
perimeter2_uc numeric,
cmp numeric,
sliver character varying(1)
);';
RAISE NOTICE 'Find all overlappings with geometrytype Polygon. ';
EXECUTE
'INSERT INTO '||schema_a||'.pgoverlap_tmp1 (
a_'||a_id||', b_'||a_id||', a_geom, b_geom, i_geom)
SELECT DISTINCT
sel.a_'||a_id||',
sel.b_'||a_id||',
sel.a_geom,
sel.b_geom,
sel.i_geom
FROM (
SELECT
a.'||a_id||' AS a_'||a_id||',
b.'||a_id||' AS b_'||a_id||',
a.the_geom AS a_geom,
b.the_geom AS b_geom,
(ST_Dump(ST_Intersection(a.the_geom,b.the_geom))).geom AS i_geom
FROM
'||schema_a||'.'||table_a||' a
LEFT JOIN
'||schema_a||'.'||table_a||' b
ON a.the_geom && b.the_geom
WHERE
a.gid < b.gid AND
ST_Intersects(a.the_geom, b.the_geom)) AS sel
WHERE
ST_GeometryType(sel.i_geom) NOT IN (
''ST_Point'',
''ST_LineString'',
''ST_MultiPoint'',
''ST_MultiLineString'',
''ST_Line'');';
EXECUTE
'CREATE INDEX pgoverlap_tmp1_gid_btree
ON '||schema_a||'.pgoverlap_tmp1
USING btree(gid);';
EXECUTE
'CREATE INDEX pgoverlap_tmp1_a_'||a_id||'_btree
ON '||schema_a||'.pgoverlap_tmp1
USING btree(a_'||a_id||');';
EXECUTE
'CREATE INDEX pgoverlap_tmp1_b_'||a_id||'_btree
ON '||schema_a||'.pgoverlap_tmp1
USING btree(b_'||a_id||');';
--Check for sliver polygons
RAISE NOTICE 'Check for sliver polygons. ';
--1. Calculation of area and perimeter:
EXECUTE
'UPDATE '||schema_a||'.pgoverlap_tmp1
SET area=ST_Area(i_geom),
perimeter2=ST_Perimeter(i_geom)^2;';
--2. Calculation of circle perimeter with area equal to area of the polygon:
EXECUTE
'UPDATE '||schema_a||'.pgoverlap_tmp1
SET perimeter2_uc=4*pi()*area;';
--3. Calculation of Compactness Index (cmp):
EXECUTE
'UPDATE '||schema_a||'.pgoverlap_tmp1
SET cmp=perimeter2_uc/perimeter2;';
--4. Decision: sliver polygon, true or false?:
EXECUTE
'UPDATE '||schema_a||'.pgoverlap_tmp1 SET sliver=''t'' WHERE area < 10;';
EXECUTE
'UPDATE '||schema_a||'.pgoverlap_tmp1 SET sliver=''t'' WHERE cmp < 0.13;';
EXECUTE
'UPDATE '||schema_a||'.pgoverlap_tmp1 SET sliver=''f'' WHERE sliver IS NULL;';
/*
EXECUTE
'UPDATE '||schema_a||'.pgoverlap_tmp1 SET sliver=''f'' WHERE area > 150;';
*/
--Create table for sliver polygons:
EXECUTE
'DROP TABLE IF EXISTS '||schema_a||'.pgoverlap_tmp2;';
EXECUTE
'CREATE TABLE '||schema_a||'.pgoverlap_tmp2 (
gid serial,
a_'||a_id||' integer,
a_geom geometry,
i_geom_union geometry,
the_geom geometry
);';
RAISE NOTICE 'Processing sliver polygons. ';
EXECUTE
'INSERT INTO '||schema_a||'.pgoverlap_tmp2
(a_'||a_id||', a_geom, i_geom_union)
SELECT a_'||a_id||', a_geom, CASE WHEN ST_Union(i_geom) IS NULL THEN ST_Collect(i_geom) ELSE ST_Union(i_geom) END
FROM '||schema_a||'.pgoverlap_tmp1
WHERE sliver=''t''
GROUP BY a_'||a_id||', a_geom;';
i:=0;
sql_1:='SELECT gid FROM '||schema_a||'.pgoverlap_tmp2 ORDER BY gid;';
FOR i IN EXECUTE sql_1 LOOP
BEGIN
EXECUTE
'UPDATE '||schema_a||'.pgoverlap_tmp2 SET the_geom=ST_Difference(a_geom, i_geom_union) WHERE gid='||i||';';
EXCEPTION WHEN internal_error THEN
EXECUTE
'UPDATE '||schema_a||'.pgoverlap_tmp2 SET the_geom=ST_Difference(ST_Buffer(a_geom,0.0), ST_Buffer(i_geom_union,0.0)) WHERE gid='||i||';';
END;
END LOOP;
EXECUTE
'CREATE INDEX pgoverlap_tmp2_a_'||a_id||'_btree
ON '||schema_a||'.pgoverlap_tmp2
USING btree(a_'||a_id||');';
--Create table with IDs of sliver polygons:
RAISE NOTICE 'Creating table with IDs of sliver polygons. ';
EXECUTE
'DROP TABLE IF EXISTS '||schema_a||'.pgoverlap_tmp3;';
EXECUTE
'CREATE TABLE '||schema_a||'.pgoverlap_tmp3 (
'||a_id||' integer);';
EXECUTE
'INSERT INTO '||schema_a||'.pgoverlap_tmp3 ('||a_id||')
SELECT DISTINCT a_'||a_id||'
FROM '||schema_a||'.pgoverlap_tmp2;';
EXECUTE
'INSERT INTO '||schema_a||'.pgoverlap_tmp3 ('||a_id||')
SELECT DISTINCT a.b_'||a_id||'
FROM
'||schema_a||'.pgoverlap_tmp1 a
LEFT JOIN
'||schema_a||'.pgoverlap_tmp2 b
ON a.b_'||a_id||'=b.a_'||a_id||'
WHERE b.a_'||a_id||' IS NULL;';
--Insert non-sliver polygons into separate table:
BEGIN
EXECUTE
'SELECT DropGeometryTable('''||schema_a||''','''||table_result2||''');';
EXCEPTION WHEN undefined_table THEN
END;
EXECUTE
'CREATE TABLE '||schema_a||'.'||table_result2||' (
gid serial PRIMARY KEY,
a_'||a_id||' integer,
b_'||a_id||' integer
);';
EXECUTE
'SELECT AddGeometryColumn('''||schema_a||''', '''||table_result2||''',''a_geom'',(SELECT DISTINCT ST_SRID(the_geom) FROM '||schema_a||'.'||table_a||'),''POLYGON'',2);';
EXECUTE
'SELECT AddGeometryColumn('''||schema_a||''', '''||table_result2||''',''b_geom'',(SELECT DISTINCT ST_SRID(the_geom) FROM '||schema_a||'.'||table_a||'),''POLYGON'',2);';
EXECUTE
'SELECT AddGeometryColumn('''||schema_a||''', '''||table_result2||''',''i_geom'',(SELECT DISTINCT ST_SRID(the_geom) FROM '||schema_a||'.'||table_a||'),''POLYGON'',2);';
RAISE NOTICE 'Inserting non-sliver polygons into separate table. ';
EXECUTE
'INSERT INTO '||schema_a||'.'||table_result2||' (
a_'||a_id||',
b_'||a_id||',
a_geom,
b_geom,
i_geom)
SELECT
a_'||a_id||',
b_'||a_id||',
a_geom,
b_geom,
i_geom
FROM '||schema_a||'.pgoverlap_tmp1
WHERE sliver=''f'';';
--create table with IDs of non-sliver polygons:
RAISE NOTICE 'creating table with IDs of non-sliver polygons. ';
EXECUTE
'DROP TABLE IF EXISTS '||schema_a||'.pgoverlap_tmp4;';
EXECUTE
'CREATE TABLE '||schema_a||'.pgoverlap_tmp4 (
'||a_id||' integer);';
EXECUTE
'INSERT INTO '||schema_a||'.pgoverlap_tmp4 ('||a_id||')
SELECT DISTINCT a_'||a_id||'
FROM '||schema_a||'.'||table_result2||';';
EXECUTE
'INSERT INTO '||schema_a||'.pgoverlap_tmp4 ('||a_id||')
SELECT DISTINCT a.b_'||a_id||'
FROM
'||schema_a||'.'||table_result2||' a
LEFT JOIN
'||schema_a||'.'||table_result2||' b
ON a.b_'||a_id||'=b.a_'||a_id||'
WHERE b.a_'||a_id||' IS NULL;';
--Insert geometries without overlappings into new table:
EXECUTE
'DROP TABLE IF EXISTS '||schema_a||'.pgoverlap_tmp5;';
EXECUTE
'CREATE TABLE '||schema_a||'.pgoverlap_tmp5 (
gid serial PRIMARY KEY,
'||a_id||' integer,
the_geom geometry
);';
EXECUTE
'INSERT INTO '||schema_a||'.pgoverlap_tmp5
('||a_id||', the_geom)
SELECT
a.'||a_id||',
a.the_geom
FROM
'||schema_a||'.'||table_a||' a
LEFT JOIN
'||schema_a||'.pgoverlap_tmp3 b
ON a.'||a_id||'=b.'||a_id||'
LEFT JOIN
'||schema_a||'.pgoverlap_tmp4 c
ON a.'||a_id||'=c.'||a_id||'
WHERE
b.'||a_id||' IS NULL AND
c.'||a_id||' IS NULL;';
--Create and merge result table:
BEGIN
EXECUTE
'SELECT DropGeometryTable('''||schema_a||''','''||table_result1||''');';
EXCEPTION WHEN undefined_table THEN
END;
EXECUTE
'CREATE TABLE '||schema_a||'.'||table_result1||' (
gid serial PRIMARY KEY,
'||a_id||' integer
);';
EXECUTE
'SELECT AddGeometryColumn('''||schema_a||''','''||table_result1||''',''the_geom'',(SELECT DISTINCT ST_SRID(the_geom) FROM '||schema_a||'.'||table_a||'),''POLYGON'',2);';
-- 11May12 PostGIS 2.0 rejects this DROP CONSTRAINT as an error -dbb
EXECUTE
'ALTER TABLE '||schema_a||'.'||table_result1||' DROP CONSTRAINT enforce_geotype_the_geom;';
--Insert Difference (sliver overlappings) into result table:
EXECUTE
'INSERT INTO '||schema_a||'.'||table_result1||' (
'||a_id||',
the_geom
)
SELECT
a_'||a_id||',
(ST_Dump(the_geom)).geom
FROM
'||schema_a||'.pgoverlap_tmp2
WHERE ST_IsEmpty(the_geom)=''f'';';
--Insert unchanged polygons (sliver overlappings) into result table:
EXECUTE
'INSERT INTO '||schema_a||'.'||table_result1||' (
'||a_id||',
the_geom
)
SELECT DISTINCT
a.b_'||a_id||',
a.b_geom
FROM
'||schema_a||'.pgoverlap_tmp1 a
LEFT JOIN
'||schema_a||'.pgoverlap_tmp2 b
ON a.b_'||a_id||'=b.a_'||a_id||'
WHERE
a.sliver=''t'' AND
b.a_'||a_id||' IS NULL;';
--Insert non-overlapping polygons into result table:
EXECUTE
'INSERT INTO '||schema_a||'.'||table_result1||' (
'||a_id||',
the_geom
)
SELECT
'||a_id||',
the_geom
FROM '||schema_a||'.pgoverlap_tmp5;';
--Update geometries already processed as sliver overlappings which also have non-sliver overlappings:
EXECUTE
'UPDATE '||schema_a||'.'||table_result2||' t SET a_geom=a.geom FROM (SELECT '||a_id||', CASE WHEN ST_Union(the_geom) IS NULL THEN ST_Collect(the_geom) ELSE ST_Union(the_geom) END AS geom FROM '||schema_a||'.'||table_result1||' GROUP BY '||a_id||') a WHERE t.a_'||a_id||'=a.'||a_id||';';
EXECUTE
'UPDATE '||schema_a||'.'||table_result2||' t SET b_geom=a.geom FROM (SELECT '||a_id||', CASE WHEN ST_Union(the_geom) IS NULL THEN ST_Collect(the_geom) ELSE ST_Union(the_geom) END AS geom FROM '||schema_a||'.'||table_result1||' GROUP BY '||a_id||') a WHERE t.b_'||a_id||'=a.'||a_id||';';
EXECUTE
'DELETE FROM '||schema_a||'.'||table_result1||' WHERE '||a_id||' IN (SELECT '||a_id||' FROM '||schema_a||'.pgoverlap_tmp4);';
--Remove tmp tables:
EXECUTE
'DROP TABLE '||schema_a||'.pgoverlap_tmp1;';
EXECUTE
'DROP TABLE '||schema_a||'.pgoverlap_tmp2;';
EXECUTE
'DROP TABLE '||schema_a||'.pgoverlap_tmp3;';
EXECUTE
'DROP TABLE '||schema_a||'.pgoverlap_tmp4;';
EXECUTE
'DROP TABLE '||schema_a||'.pgoverlap_tmp5;';
END;
$BODY$
LANGUAGE 'plpgsql' VOLATILE
COST 100;
ALTER FUNCTION public._pgoverlap(varchar, varchar, varchar, varchar, varchar) OWNER TO postgres;
