[postgis-users] st_difference() problem?

[Chris Hermansen]

Ladies and gentlemen;

Nice to be back.

I am wrestling with st_difference() and it is getting the better of me.  I
am hoping someone who has seen this before can help.

I have two tables that are buffer polygons, built by PostGIS, around some
linear features (trails and traplines).

I am seeking to combine these two tables to produce a single set of polygons
that shows areas that are inside the trail buffers, inside the trapline
buffers, both, or neither.

The pattern I am following is

st_difference(btrail, btrapline)
union
st_intersection(btrail, btrapline)
union
st_difference(btrapline, btrail)

The btrail table has 9 rows, all valid POLYGONs; the btrapline table has 18
rows, all valid POLYGONs.  Nothing self-overlaps.  All have SRID of 3005.

There is a place in the southwest of the area where one trapline buffer and
one trail buffer cross over each other twice.  No other features come close.

If I copy just those two polygons into two separate test tables and use my
difference/intersection/difference operation, I get the result I expect - 10
polygons, none overlapping, where some are trapline-only, some are
trail-only, and some are both trail and trapline.  A visual inspection of
the result in QGIS "makes sense".

If I perform the same analysis on the full data set, I get a result that
seems wrong - looking at the portion of the data related to those two
features, instead of 10 polygons, I see 35 polygons.  Many of these "result"
polygons overlap each other.  Visually (in QGIS again), I see what appears
to be the two original input polygons covering a bunch of other polygons.

I can of course post the data; I also have a few screen shots from QGIS
saved as jpg files.

This is from PostgreSQL 8.4.8 running on Ubuntu 10.10.  The output from
postgis_full_version() is

POSTGIS="1.5.1" GEOS="3.2.0-CAPI-1.6.0" PROJ="Rel. 4.7.1, 23 September 2009"
LIBXML="2.7.6" USE_STATS

This is the code that produces the multiple overlapping result polygons with
the full data set and seemingly normal results with the two test buffers.

-- create the diff/int/diff on the full table

drop table utrailtrapline;

create table utrailtrapline (
    gid serial,
    istrail boolean,
    istrapline boolean);

SELECT AddGeometryColumn('','utrailtrapline','the_geom','3005','POLYGON',2);

insert into utrailtrapline (istrail, istrapline, the_geom)
select
    istrail,
    FALSE as istrapline,
    (st_dump(st_difference(btrail.the_geom, btrapline.the_geom))).geom as
the_geom
from btrail, btrapline;

insert into utrailtrapline (istrail, istrapline, the_geom)
select
    istrail,
    istrapline,
    (st_dump(st_intersection(btrail.the_geom, btrapline.the_geom))).geom as
the_geom
from btrail, btrapline
where st_intersects(btrail.the_geom, btrapline.the_geom);

insert into utrailtrapline (istrail, istrapline, the_geom)
select
    FALSE as istrail,
    istrapline,
    (st_dump(st_difference(btrapline.the_geom, btrail.the_geom))).geom as
the_geom
from btrail, btrapline;

Thanks in advance!

-- 
Chris Hermansen
TECO Natural Resource Group Limited
301 · 958 West 8th Avenue
Vancouver BC CANADA · V5Z 1E5
Tel +1.604.714.2878 · Cel +1.778.840.4625
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20110530/7f6031f9/attachment.html>