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

[Martin Davis]

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
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>

-- 
Martin Davis
Senior Technical Architect
Refractions Research, Inc.
(250) 383-3022