A project for CE 5383 (GIS Analyses) at UC Denver involved creating a website to collect information about activities on a piece of land. Users would log in to the website, look at a map and draw some polygons representing the locations of different activities. A custom PostGIS SQL function had to be created to calculate the intersections of these polygons and to report the areas of the intersections. The purpose of this was to find use conflicts and use intensities. The custom function used PostGIS’s built-in intersection function, st_intersect. The st_intersect function only works on two geometries simultaneously, so it had to be used recursively to solve the problem.
The following movie is a demonstration of the web application using the custom PL/SQL funtion written to handle polygon intersections recursively.
The current PostGIS (1.2.1) intersection function (see PostGIS Reference 6.1.3) only handles 2 input geometries. The custom function uses the intersection function recursively over any number of polygons. The basic algorithm is:
- Place all source polygons in a new temporary table
- Pair up the polygons in the temporary table and see which ones intersect each other
- Take the first intersection and add children to the temporary table and remove parents
- Repeat 2 – 3 until there are no more intersections in the temporary table
The following slides show this algorithm applied to a simple example. Press the Next and Previous buttons to move through the slide
The function works OK with a few simple polygons, but splinters are a problem. If I had to do this again I'd convert the polygons to rasters and use matrix algebra.
And now for the actual code, set up some preliminary tables first. Testpolygons will hold the source polygons. Insert 3 polygons into this table. Then create table to hold the results, testpolygons_temp
create table testpolygons(
gid integer,
the_geom geometry,
CONSTRAINT pk PRIMARY KEY(gid)
)
WITHOUT OIDS;
insert into testpolygons (gid,the_geom) values(1, GeomFromText('POLYGON((-104.90 39.89,-105.04 39.83,-105.01 39.77,-104.85 39.80,-104.90 39.89))') ) ;
insert into testpolygons (gid,the_geom) values(2, GeomFromText('POLYGON((-105.11 39.76,-105.04 39.82,-104.97 39.80,-104.96 39.73,-105.06 39.71,-105.11 39.76))') );
insert into testpolygons (gid,the_geom) values(3, GeomFromText('POLYGON((-105.09 39.79,-105.03 39.76,-105.00 39.82,-105.04 39.87,-105.10 39.84,-105.10 39.84,-105.09 39.79))') );
CREATE TABLE testpolygons_temp
(
gid integer,
the_geom geometry,
count integer
CONSTRAINT pk2 PRIMARY KEY(gid)
)
WITH OIDS;
The first function written is a_intersect_recursive(). This will be called like so: select a_intersect_recursive(); I like to put the 'a_' in the function name so it shows up at the top of the function list.
CREATE OR REPLACE FUNCTION a_intersect_recursive () RETURNS boolean AS
$BODY$
DECLARE
polygon_id int;
k int;
the_polygons int[];
first_geom geometry;
second_geom geometry;
dynamic_query varchar;
first_gid int;
second_gid int;
int_avail boolean;
BEGIN
--create a temporary table that will initially have all the original polygons but will gradually be replaced with their intersection fragments
DROP TABLE testpolygons_temp;
CREATE TABLE testpolygons_temp AS (SELECT gid, the_geom , 1 as count FROM testpolygons);
EXECUTE 'ALTER TABLE testpolygons_temp ADD PRIMARY KEY (gid)';
--go through the array and create 1:1 pairs between the gids. see if any of those polygons intersect each other
--if you had polygons with ids of 1,2,3,4,5 this would look at 12, 13, 14, 15, 23, 24,25, 34, 35, and 45
--take the first intersection and explode those polygons into their fragments . erase the originals
--find the first intersection
int_avail = true;
WHILE (int_avail) LOOP
--add the gid of every polygon in testpolygons into an array
the_polygons = ARRAY[0];
int_avail = false;
EXECUTE 'SELECT ARRAY(SELECT gid FROM testpolygons_temp )' INTO the_polygons;
RAISE NOTICE 'array: % ', the_polygons;
FOR j IN 1..array_upper(the_polygons, 1)LOOP
FOR k IN 1..array_upper(the_polygons, 1)LOOP
IF (j != k) THEN
first_gid = the_polygons[j];
second_gid = the_polygons[k];
EXECUTE 'SELECT the_geom from testpolygons_temp where gid='||first_gid INTO first_geom;
EXECUTE 'SELECT the_geom from testpolygons_temp where gid='||second_gid INTO second_geom;
--exit the inner loop if there is an intersection
IF intersects(first_geom,second_geom) THEN
--keep small slivers out by keeping the area of the intersection be more than {something more than 0}
IF area(intersection(first_geom,second_geom))>0.000000001 THEN
IF NumGeometries(first_geom) is null THEN
IF NumGeometries(second_geom) is null THEN
IF isvalid(first_geom) THEN
IF isvalid(second_geom) THEN
int_avail = true;
EXECUTE 'SELECT a_intersect_and_save((SELECT the_geom FROM testpolygons_temp WHERE gid='||first_gid||'), (select the_geom from testpolygons_temp where gid='||second_gid||'), (SELECT count FROM testpolygons_temp WHERE gid='||first_gid||'), (SELECT count FROM testpolygons_temp WHERE gid='||second_gid||') )';
EXECUTE 'DELETE FROM testpolygons_temp WHERE gid ='||first_gid;
EXECUTE 'DELETE FROM testpolygons_temp WHERE gid ='||second_gid;
END IF;
END IF;
END IF;
END IF;
END IF;
END IF;
END IF;
END LOOP;
END LOOP;
END LOOP;
RETURN TRUE;
END;
$BODY$
LANGUAGE 'plpgsql' VOLATILE;
The second function is a_intersect_and_save. This function is called over and over by the first function for each pair of polygons. The input geometries are buffered with 0 units. This seems to cause less errors.
CREATE OR REPLACE FUNCTION a_intersect_and_save(poly1 geometry, poly2 geometry, poly1count int, poly2count int)
RETURNS VOID AS
$BODY$
DECLARE
poly1small geometry;
poly2small geometry;
the_intersection geometry;
dynamic_query1 varchar;
dynamic_query2 varchar;
dynamic_query3 varchar;
sum_of_counts int;
BEGIN
-- find the intersection between the two polygons. this will become the first child.
the_intersection := intersection(buffer(poly1,0), buffer(poly2,0));
--RAISE NOTICE 'Intersection area: % ', area(the_intersection);
IF area(the_intersection) = 0 THEN
RETURN;
END IF;
--RAISE NOTICE 'poly and intersectoin % ', overlaps(poly1, the_intersection);
--find the difference between polygon 1 and polygon 2. this will become a new child polygon
poly1small := difference(buffer(poly1,0), buffer(poly2,0));
--RAISE NOTICE 'poly1 minus intersection area: % ', area(poly1small);
--find the difference between polygon 2 and polygon 1. this will become a new child polygon
poly2small := difference(buffer(poly2,0), buffer(poly1,0));
sum_of_counts = poly1count + poly2count;
dynamic_query1 = E'INSERT INTO testpolygons_temp values ((select max(gid) from testpolygons_temp)+1, \'';
dynamic_query1 = dynamic_query1 || the_intersection ;
dynamic_query1 = dynamic_query1 || E'\',';
dynamic_query1 = dynamic_query1 || sum_of_counts || ')';
EXECUTE dynamic_query1;
dynamic_query2 = E'INSERT INTO testpolygons_temp values ((select max(gid) from testpolygons_temp)+1, \'';
dynamic_query2 = dynamic_query2 || poly1small ;
dynamic_query2 = dynamic_query2 || E'\',';
dynamic_query2 = dynamic_query2 || poly1count || ')';
EXECUTE dynamic_query2;
dynamic_query3 = E'INSERT INTO testpolygons_temp values ((select max(gid) from testpolygons_temp)+1, \'';
dynamic_query3 = dynamic_query3 || poly2small ;
dynamic_query3 = dynamic_query3 || E'\',';
dynamic_query3 = dynamic_query3 || poly2count || ')';
EXECUTE dynamic_query3;
RETURN;
END;
$BODY$
LANGUAGE 'plpgsql' VOLATILE;
And that's it! I like to use QuantumGIS as a viewer connected to testpolygons and testpolygons_temp.