[postgis-users] Chained Intersections?

[Curtis W. Ruck]

This is a requirement given to me.  (I had no say in it, but it would make my application ultra cool)

I need to have the ability to output from the database all possible intersections over a set of polygons.  I wrote a pl/pgsql function Intersection(GeometryCollection) that iterates over the geometry collection and returns the intersection of all given polygons.  The problem is that if all of the polygons don't share an area then it returns nothing.  I need to write a pl/pgsql function that takes a GeometryCollection and outputs a set of polygons and integers.  Each polygon would have an integer showing how many polygons intersected to create that polygon.

For example, if you look at http://orbyn.mine.nu/intersections/4polys.png (to save the mailing list the attachment) the input polygons (red, green, blue) should return (yellow, cyan, purple, white) from the function where yellow, cyan, purple have a value of 2 and white has a value of 3.

I wrote a pl/pgsql function that i thought would do this, far below labeled Intersection_All.  The problem with this is that it needs to be called from a FROM clause because it returns a SETOF records, but within a FROM clause i can't do an aggregate method collect(the_geom).  (I also tried a subquery in the FROM clause).

Does anyone have any suggestions on how to warp Postgres/PostGIS into doing this?

Curtis W. Ruck
<unnamed large entity>



CREATE OR REPLACE FUNCTION intersection_all(in_polys geometry, OUT intersection_count int, OUT the_geom geometry)
  RETURNS SETOF record AS
$BODY$
DECLARE
	total_count int := 0;
	current_level_count int := 0;
	first_id_of_previous_level int := 0;
	last_id_of_previous_level int:= 0;
	num_intersections int ARRAY[0];
	geometries geometry ARRAY[0];
	temp_geom geometry;
BEGIN
	FOR i IN 1..NumGeometries(in_polys) LOOP
		FOR j IN i+1..NumGeometries(in_polys) LOOP
			IF Intersects(GeometryN(in_polys,i),GeometryN(in_polys,j)) THEN
				temp_geom := Intersection(GeometryN(in_polys,i),GeometryN(in_polys,j));
				geometries := geometries || temp_geom;
				num_intersections := num_intersections || 1;
				last_id_of_previous_level := last_id_of_previous_level+1;
				total_count := total_count +1;
			END IF;
		END LOOP;
	END LOOP;
 
	FOR level IN 2..100 LOOP
		FOR i IN first_id_of_previous_level..last_id_of_previous_level LOOP
			FOR j IN i+1..last_id_of_previous_level LOOP
				IF Intersects(geometries[i],geometries[j]) THEN
					temp_geom := Intersection(geometries[i],geometries[j]);
					current_level_count := current_level_count+1;
					geometries := geometries || temp_geom;
					num_intersections := num_intersections || level;
					total_count := total_count +1;
				END IF;
			END LOOP;
		END LOOP;
		IF current_level_count > 0 THEN
			first_id_of_previous_level := last_id_of_previous_level+1;
			last_id_of_previous_level := first_id_of_previous_level+current_level_count;
			current_level_count := 0;
		ELSE 
			EXIT;
		END IF;
	END LOOP;
 
	FOR i IN 1..total_count LOOP
		the_geom = geometries[i];
		intersection_count = num_intersections[i];
		RETURN NEXT;
	END LOOP;
	RETURN;
END;
$BODY$
  LANGUAGE 'plpgsql' IMMUTABLE;



-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20061009/39b47680/attachment.html>