[postgis-users] Re: computing unions of intersecting polygons

[Reid Priedhorsky]

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));