[postgis-users] PostGIS analogues to ArcInfo geoprocessing functions

[Martin Davis]

The big problem with this in PostGIS land is that Union (or let's call 
it Overlay, to distinguish it from the existing Union) is a *global* 
function.  I.e. it operates over *all* polygons in a table (or set of 
tables).  This is in contrast to most existing PostGIS functions, which 
operate over a single or a pair of Geometries

Your example is a bit simplistic - typically people want to overlay 
1000's of polygons at a time.  There's two ramification of this. 

(a) This can't be done in a piecewise, polygon-by-polygon fashion (or at 
least, not efficiently, correctly, or robustly). 
(b) This can't be done in-memory.  It requires some external structures 
(in PostGIS, presumably tables, although possibly some other file-based 
structure would work as well or better).

ArcINFO must have some long-standing rocket science to compute 
overlays.  This is majorly helped by the fact that its native data 
structure is an arc-node topology, which is a necessary precursor to 
computing an overlay efficiently - and which is itself expensive and 
difficult to compute in a non-memory bound fashion.

The basic approach to computing polygon overlays has been 
well-understood for a long time (although this does not imply 
well-documented!).  The implementation however is quite tricky, 
especially if performance and robustness is required.  There is a 
notable lack of open-source existing implementations - which should be 
an indication of just how difficult this really is.

Not that it wouldn't be fun to try... but this is a multi-month effort 
for some clever developers familiar with the area.



Chris Hermansen wrote:
> Following on from Paul's note about a function that would make the 
> output of st_intersection() more immediately useful, I have a 
> favourite whine to mention.
>
> Some geoprocessing activities that are easy to do in ArcInfo are 
> more-or-less as easy to do in PostGIS.  One, though, seems to me to be 
> very difficult or slow or both, and that is the ArcInfo UNION function.
>
> A simple example:  Let's say I have polygon theme "a" as (0 0, 2 0, 2 
> 1, 0 1, 0 0) and "b" as (1 0, 3 0, 3 1, 1 1, 1 0).  Then
>
> UNION a b c
>
> will produce polygon theme "c" with three polygons: (0 0, 1 0, 1 1, 0 
> 1, 0 0), (1 0, 2 0, 2 1, 1 1, 1 0), and (2 0, 3 0, 3 1, 2 1, 2 0).
>
> In other words, the three polygons in "c" bound: the points unique to 
> the polygon in "a"; the points common to the polygons in "a" and "b"; 
> and the points unique to the polygon in "b".
>
> This seems exceedingly difficult in PostGIS; for example, see the 
> article at 
> http://www.postgis.org/support/wiki/index.php?ExamplesOverlayTables
>
> I can vouch from personal experience, at least following the approach 
> recommended at the above link by Kevin, that this is not a 
> particularly speed operation either.
>
> And yet, so easy to do in ArcInfo.  And, a pretty useful and basic 
> kind of geoprocessing function.  In my mind, it is a kind of spatial 
> outer join, where output polygons have either the attributes from "a", 
> the attributes from "b", or both.
>
> Another way to think of this is to use the polygons from "b" to make 
> holes in "a", putting the results in "c1".  Then use the polygons from 
> "a" to make holes in "b", putting the results in "c2".  Finally, 
> intersect "a" and "b" and put that in "c3".  Then append c1, c2, and 
> c3 together.
>
> Please tell me I'm not the only one "out there" who's stuck on this 
> problem; or, if I am, please tell me that function 
> "st_something_i_havent_tried_yet()" that gives me the result I'm 
> looking for - perhaps wrapped in "st_just_gimme_the_polygons()".
>

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