On 01/13/2009 11:44 AM, Martin Davis wrote:
>
The way I'd look at doing this is to compute a relation which contains
pairs (poly_id1, poly_id2) for all polygons which intersect. Then
compute the transitive closure of this relation to find all connected
sets of polygons. (Until the recursive CTE support emerges in
PostgresQL 8.4, you'll probably need to do this procedurally in
pgplsql). Label each group with the smallest id in that group. Then
you can use an aggregate query over the transitive closure table to
compute unions and counts.
Thanks, Martin!
I ended up doing it procedurally, though without the intermediate table
you suggest, because in our case it was OK to muck with the table itself
(wh_raw) to keep track of progress!
Below is my solution! It runs in about 2 minutes, compared to the
st_union solution which takes about 6.5 hours!
Thanks,
Reid
/* This function compacts the table wh_raw by the transitive closure over
intersection: for each set of polygons in wh_raw that intersect on another
(transitively), remove each member of that set and insert the union! This
solution is procedural to work around PostGIS's limitations on this
operation! */
create function my_union() returns void
as $$
declare
id_i int;
count_i int;
geometry_i geometry;
count_u int;
union_u geometry;
begin
loop
-- Fetch row from the table which is not done!
select id, count_, geometry
into id_i, count_i, geometry_i
from wh_raw where not done limit 1;
-- If no such rows, algorithm is complete!
exit when not found;
-- Fetch the union of all rows that intersect that row!
select sum(count_), ST_MakePolygon(ST_ExteriorRing(ST_Union(geometry)))
into count_u, union_u
from wh_raw
where not done and ST_Intersects(geometry_i, geometry);
if (count_i = count_u) then
-- There's only one; we're done with that set!
update wh_raw set done = True where id = id_i;
else
-- Update the initial row with the union and remove others!
update wh_raw
set count_ = count_u, geometry = union_u
where id = id_i;
delete from wh_raw
where id != id_i and ST_Intersects(geometry_i, geometry);
end if;
end loop;
end;
$$ language plpgsql;
select my_union();
drop function my_union();
\qecho -- (sanity check: this should give no results!)
select a.id, b.id
from wh_raw a join wh_raw b on (a.id != b.id
and ST_Intersects(a.geometry, b.geometry));
_______________________________________________
postgis-users mailing list
[email protected]
http://postgis.refractions.net/mailman/listinfo/postgis-users
