[postgis-users] Re: computing unions of intersecting polygons
Reid Priedhorsky
reid at umn.edu
Wed Jan 14 15:24:54 PST 2009
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));
More information about the postgis-users
mailing list