[geos-devel] Small gaps between results of GEOSDifference?

[Sandro Santilli]

On Fri, Feb 15, 2013 at 02:40:16PM +0100, Marco Hugentobler wrote:
> >If you want to further
> >debug your case you may try to build geos yourself with uncommenting
> >the GEOS_DEBUG_BINARYOP define in include/geos/geom/BinaryOp.h and
> >see what happens to each operation.
> 
> Ok, the output from the program is now:
> 
> Trying with original input.
> Trying with original input.
> Trying with original input.
> POLYGON ...
> POLYGON ...
> 
> Does 'Trying with original input' mean there are no robustness issues?

It means no exception has been thrown.
So the problem is NOT with the heuristics in BinaryOp.

After taking a look at the geometries I now know what's going on.
The two very simple (almost rectangle) polygons overlap so that
their boundary intersect at exactly 2 points.

Those 2 points are NOT vertices in the input geometries, so when
constructing a geometry resulting from an overlay (be it union
or difference) a vertex is added to the boundary producing a
slight modification of the so-splitted segment.
This is expected and unavoidable with finite precision.

In the A + B case all geometries take part to the build of a single
topology, so the intersection points become vertexes of the single
output geometry.

In the A + ( B - A ) instead there are two different topologies
being built. First topology (B-A) is built the same way as (B+A)
would be built, produces the same vertices and gives you back
a polygon with those two additional vertices that are NOT NECESSARELY
found on the segment they did split.

The second topology building phase (A + (B-A)) finds one of those
two vertices is NOT on the line, and so does not generate noding at
that specific point.

I see two possible solutions to this problem:

 1. Use arbitrary precision
 2. Use controlled fixed precision

GEOS doesn't implement 1, and I don't think that's generally a good solution.
Fixed precision _is_ implemented BUT not exposed to the C-API yet.

What you could do, as a workaround, is explicitly snapping the geometries
each other to a tolerance of choice. That way geometry A would snap to
the new node - now being a vertex of (B-A) - and topology builder would 
node them consequently.

Exposing fixed precision to the C-API is on the TODO list for 3.1 so if
you want to help with that you're very welcome. It will take exposing
a new PrecisionModel abstract class and provide means to attach them to
GEOSGeometry objects.

See also: http://trac.osgeo.org/postgis/wiki/ToleranceDiscussion

--strk; 

 http://www.cartodb.com - Map, analyze and build applications with your data

                                       ~~ http://strk.keybit.net