[postgis-users] PostGIS analogues to ArcInfo geoprocessing functions
Martin Davis
mbdavis at refractions.net
Thu Jan 29 08:42:31 PST 2009
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
More information about the postgis-users
mailing list