[Gdal-dev] Re: ogr geometry intersection

[Todd Jellett]

Hi Jack,

You have to be very careful when mixing "tests" (ie the binary 
predicates disjoint, contains, ...) with "constructions" (the set 
theoreticals intersection, union, ...). This will work most of the time 
but not always. This is because gdal's geometry backend is geos and geos 
uses infinite precision for the tests and a finite precision for 
constructions. This means that there will be cases where disjoint() and 
friends returns true but if you go ahead and do the intersection, you 
will find that the intersection is finite (non-NULL) and valid.

This discontinuity between the tests and constructions makes it very 
difficult to implement mission critical GIS software (you can read 
engineering for $$$ from this) with geos. This is because whenever you 
have a construct like:

if ( g1.disjoint( g2 ) == false )  // or any of the other tests
{
   g3 = g1.intersection( g2 );
    ...
}

there is a possibility that the program flow will fail the test (the if 
clause) and not do an intersection when one should be done and anytime 
the program flow fails to execute a block of code because the the test 
in the if clause fails, this is a bug (IMHO).

The cases that generate this kind of failure come up in (our) GIS 
software more frequently than you would think. We had to re-write our 
own (inhouse) geometries (with tests equal to constructions) several 
years ago (2000).

Todd

Jack Riley wrote:
> It appears that the FWTools v1.0.9-and-earlier TopologyException is 
> limited to situations where one of the two polynomials being tested 
> for Intersection() is completely within the other.  In other words, 
> the Intersection() test in "old" versions of FWTools/GDAL works--and 
> works "fast"--if the two polynomials actually intersect.  So, my 
> workaround is to use FWTools v1.0.9 (or earlier) and condition any 
> call to Intersection() with Disjoint() or Contains() test(s).  The 
> question remaining is:  can FWTools v1.1.0+/GDAL revert back to the 
> use of the "old" and "fast" Intersection() code, adding in a test of 
> Disjoint()/Contains() to circumvent the TopologyException?
>
> On 3/21/07, *Jack Riley* <jack.l.riley at gmail.com 
> <mailto:jack.l.riley at gmail.com>> wrote:
>
>     Hello,
>
>     I recently upgraded my FWTools (Win32) installation from v1.0.5 to
>     v1.2.2.  Using v1.0.5, ogr geometry intersections would
>     occasionally raise a TopologyException in my Python application; e.g.:
>
>     Original exception: TopologyException: Area of intersection result
>     is bigger then minimum area between input geometries
>     Trying with Common bits remover.
>     CBR: TopologyException: Area of intersection result is bigger then
>     minimum area between input geometries
>     Trying with precision 25
>     Reduced with precision (25): TopologyException: no outgoing
>     dirEdge found - 71.6 40.88 40.8829
>     Trying with precision 24
>     Reduced with precision (24): TopologyException: no outgoing
>     dirEdge found -71.5833 40.875 40.8829
>     Trying with precision 23
>
>     After updating to FWTools v1.2.2 (gdal 1.4.0+svn to Feb. 22,
>     2007), I no longer see TopologyExceptions from polygyon
>     intersections in my Python application -- but the performance has
>     gone way down.  Evaluation of FWTools versions between 1.0.5 and
>     1.2.2 reveals that the TopologyException fix + slow-down was
>     introduced in v1.1.0:
>
>     P1.Intersection(P2), where P1=rectangle and P2=1860-point polygon
>     v1.0.9:  ~0.01 seconds
>     v1.1.0:  ~10.00 seconds (1000x slower)
>
>     See attached sample program.
>
>     Is this a bug in GDAL, FWTools (Python wrapper), or an underlying
>     library (GEOS)?
>     Any ideas re: the slow-down and possible fix are appreciated.
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> Gdal-dev mailing list
> Gdal-dev at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/gdal-dev