Beauty, eh!
Reid Priedhorsky wrote:
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
--
Martin Davis
Senior Technical Architect
Refractions Research, Inc.
(250) 383-3022
_______________________________________________
postgis-users mailing list
[email protected]
http://postgis.refractions.net/mailman/listinfo/postgis-users
