Skip to content

gIntersection returns error "TopologyException: no outgoing dirEdge found at"

15 messages · Jean-Luc Dupouey, Roger Bivand, Edzer Pebesma +1 more

Original
Jean-Luc Dupouey
# 
Why gIntersection returns the ""TopologyException: no outgoing dirEdge 
found at" error in the following very simple case:

 > #construct a spatial layer with two adjacent squares, each of 10 x 10 
units:
 >
 > Pol1=rbind(c(0,0),c(0,10),c(10,10),c(10,0),c(0,0))
 > Pol2=rbind(c(0,0),c(10,0),c(10,-10),c(0,-10),c(0,0))
 >
 > Pols1=Polygons(list(Polygon(Pol1)),"Pols1")
 > Pols2=Polygons(list(Polygon(Pol2)),"Pols2")
 > MyLay=SpatialPolygons(list(Pols1,Pols2))
 >
 > #construct the same layer, but lagged by 0.5 in x and y directions
 >
 > Pol1l=Pol1+0.5
 > Pol2l=Pol2+0.5
 >
 > Pols1l=Polygons(list(Polygon(Pol1l)),"Pols1l")
 > Pols2l=Polygons(list(Polygon(Pol2l)),"Pols2l")
 > MyLayl=SpatialPolygons(list(Pols1l,Pols2l))
 >
 > #view the resulting spatial layers:
 >
 > plot(MyLay)
 > plot(MyLayl,add=TRUE)
 >
 > #try to intersect:
 >
 > inter=gIntersection(MyLay,MyLayl)
Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, 
"rgeos_intersection") :
   TopologyException: no outgoing dirEdge found at 0.5 0

It works with the option byid=TRUE.

But my question is why it does not work without? Is this behaviour 
predictable?

I went through some previous related posts to the list, but could not 
find the answer.

Thanks,
Roger Bivand
#
On Wed, 23 Sep 2015, Jean-Luc Dupouey wrote:

            

      
      
      
    

        
The error message is coming from GEOS - you are welcome to investigate 
further. If you use gUnaryUnion() on the objects first, there is no error:

library(rgeos)
Pol1=rbind(c(0,0),c(0,10),c(10,10),c(10,0),c(0,0))
Pol2=rbind(c(0,0),c(10,0),c(10,-10),c(0,-10),c(0,0))
library(sp)
Pols1=Polygons(list(Polygon(Pol1)),"Pols1")
Pols2=Polygons(list(Polygon(Pol2)),"Pols2")
MyLay=SpatialPolygons(list(Pols1,Pols2))
Pol1l=Pol1+0.5
Pol2l=Pol2+0.5
Pols1l=Polygons(list(Polygon(Pol1l)),"Pols1l")
Pols2l=Polygons(list(Polygon(Pol2l)),"Pols2l")
MyLayl=SpatialPolygons(list(Pols1l,Pols2l))
plot(MyLay)
plot(MyLayl,add=TRUE)
points(x=0.5, y=0)
inter=gIntersection(gUnaryUnion(MyLay), gUnaryUnion(MyLayl))
plot(inter, add=TRUE, border="green")

It is a GEOS issue, not rgeos.

            

      
      
      
    

        
This is unknown. I'll try again on GEOS 3.5.0 and see if it is still 
there: yes, unfortunately it is still there.

Roger

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    
1 day later

    

      
      

        


  

        

      Jean-Luc Dupouey
    

    

      
      #
    

  


  

      
Thanks for your answer. Unfortunately, I have no time to spend debugging 
the GEOS code. Just making some tests under R, I found an even worse 
situation, where the program seems to work, but gives a completely false 
result:

library(sp)
library(rgeos)

#build a first polygon, MyLay, composed of two adjacent squares, size 1x1:

Pol1=rbind(c(0,0),c(0,1),c(1,1),c(1,0),c(0,0))
Pol2=rbind(c(0,0),c(1,0),c(1,-1),c(0,-1),c(0,0))

Pols1=Polygons(list(Polygon(Pol1)),"Pols1")
Pols2=Polygons(list(Polygon(Pol2)),"Pols2")
MyLay=SpatialPolygons(list(Pols1,Pols2))

#build a second polygon, MyLayl, composed of the same two squares lagged 
by 0.1 in x and y directions:

Pol1l=Pol1+0.1
Pol2l=Pol2+0.1

Pols1l=Polygons(list(Polygon(Pol1l)),"Pols1l")
Pols2l=Polygons(list(Polygon(Pol2l)),"Pols2l")
MyLayl=SpatialPolygons(list(Pols1l,Pols2l))

#view the resulting spatial layers:

plot(MyLay)
plot(MyLayl,add=TRUE)

#"successful" intersection:

inter=gIntersection(MyLay,MyLayl)

#but false result: the intersection gives the same contour as MyLay!

plot(MyLay)
plot(MyLayl,add=TRUE)
plot(inter,col="red",add=TRUE)

I use R version 3.2.2, rgeos_0.3-12 and sp_1.2-0.

If such an error occurs when crossing larger spatial layers, it will 
remain undetectable.

If I am not wrong, it implies that this software is not very reliable. 
Especially when considering that byid=FALSE is the default option.

gIntersection gives a correct result using the option byid=TRUE. The 
intersect function of the raster package also gives a correct 
intersection. Do you think this latter function could be a better 
option, in general?

Thanks in advance,

Jean-Luc Dupouey

INRA
Forest Ecology & Ecophysiology Unit
F-54280 Champenoux
France
mail: dupouey at nancy.inra.fr

Le 23/09/2015 19:39, Roger Bivand a ?crit :

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Fri, 25 Sep 2015, Jean-Luc Dupouey wrote:

        

            

      
      
      
    

        
Well, none of us get paid around here for support or maintenance, so you 
are the guy asking the question, you are the one interested in the answer, 
you should find the time to "pay back" what other volunteers have created.

            

      
      
      
    

        
You did not follow my advice to put the objects into gUnaryUnion() first, 
to remove the problem. It isn't obvious where this problem is coming from, 
but most likely requires debugging in GEOS itself.

            

      
      
      
    

        
This was not what you said you wanted, which was the (single) object for 
which MyLay intersected MyLayl.

            

      
      
      
    

        
No, not at all. If you inspect it, you'll see that it simply tries to 
"help" users by trying to guess which "data" slot elements might be copied 
across. It also calls gIntersect() in dense mode, which will choke on 
intersections between objects with many geometries. It runs 
gIntersection(..., byid=TRUE), so adds little to just doing that.

See:

library(rgdal)
getMethod("intersect", c("SpatialPolygons", "SpatialPolygons"))

Hope this clarifies,

Roger

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Fri, 25 Sep 2015, Roger Bivand wrote:

        

            

      
      
      
    

        
...

            

      
      
      
    

        
Changing Pol1 and Pol2 by * 5, and the increment below to 0.5 gives the 
earlier:

            

      
      
      
    

        
Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, 
"rgeos_intersection") :
   TopologyException: no outgoing dirEdge found at 0.5 0

so this is another rendering of the original edge case in GEOS in this 
thread. The practical resolution is to run gUnaryUnion() on the objects 
when byid=FALSE and the required output is a single (possibly multipart) 
geometry containing the whole intersection.

Should we consider using gUnaryUnion() by default if byid for the object 
is FALSE, but permit users to choose not to do this?

Roger

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Edzer Pebesma
    

    

      
      #
    

  


  

      
On 25/09/15 14:36, Roger Bivand wrote:

            

      
      
      
    

        
Has someone tried to reproduce this with GEOS, without R, e.g. by a
PostGIS query?

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Fri, 25 Sep 2015, Edzer Pebesma wrote:

        

            

      
      
      
    

        
Do you have a running PostGIS instance - I don't. I did try v.overlay in 
GRASS, but it isn't GEOS-based, and returns the correct answer for 
byid=c(TRUE, TRUE) so like rgeos in the same setting. I'm re-installing 
GEOS with Python, but someone else should check (my GEOS 3.5.0).

Roger

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Vijay Lulla
    

    

      
      #
    

  


  

      
Seems to be working in PostGIS!  Below is from my PostGIS session:

gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 10, 10 10, 10 0, 0
0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
 st_intersects
---------------
 t
(1 row)


gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 10, 10 10, 10 0, 0
0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
                                  st_intersection
------------------------------------------------------------------------------------
 0102000000020000000000000000002440000000000000000000000000000000000000000000000000
(1 row)

gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 1, 1 1, 1 0, 0
0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
                                  st_intersection
------------------------------------------------------------------------------------
 010200000002000000000000000000F03F000000000000000000000000000000000000000000000000
(1 row)


gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 1, 1 1, 1 0, 0
0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
 st_intersects
---------------
 t
(1 row)

gisdb=# select postgis_version();
            postgis_version
---------------------------------------
 2.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1
(1 row)


gisdb=# select postgis_full_version();

  postgis_full_vers
ion
------------------------------------------------------------------------------------------
----------------------------------------------------------------------------
 POSTGIS="2.1.7 r13414" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel.
4.8.0, 6 March 2012" GDAL
="GDAL 1.11.1, released 2014/09/24" LIBXML="2.7.8" LIBJSON="UNKNOWN" RASTER
(1 row)

        
On Fri, Sep 25, 2015 at 1:02 PM, Roger Bivand <Roger.Bivand at nhh.no> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Fri, 25 Sep 2015, Roger Bivand wrote:

        

            

      
      
      
    

        
With Shapely, I see (values output from rgeos::writeWKT):

from shapely.wkt import dumps, loads

            

      
      
      
    

        
0.0000000000000000, 0.0000000000000000 1.0000000000000000, 
1.0000000000000000 1.0000000000000000, 1.0000000000000000 
0.0000000000000000, 0.0000000000000000 0.0000000000000000)), POLYGON 
((0.0000000000000000 0.0000000000000000, 1.0000000000000000 
0.0000000000000000, 1.0000000000000000 -1.0000000000000000, 
0.0000000000000000 -1.0000000000000000, 0.0000000000000000 
0.0000000000000000)))")

            

      
      
      
    

        
0.1000000000000000, 0.1000000000000000 1.1000000000000001, 
1.1000000000000001 1.1000000000000001, 1.1000000000000001 
0.1000000000000000, 0.1000000000000000 0.1000000000000000)), POLYGON 
((0.1000000000000000 0.1000000000000000, 1.1000000000000001 
0.1000000000000000, 1.1000000000000001 -0.9000000000000000, 
0.1000000000000000 -0.9000000000000000, 0.1000000000000000 
0.1000000000000000)))")

            

      
      
      
    

        
'POLYGON ((0.0000000000000000 0.0000000000000000, 0.0000000000000000 
1.0000000000000000, 1.0000000000000000 1.0000000000000000, 
1.0000000000000000 0.0000000000000000, 1.0000000000000000 
-1.0000000000000000, 0.0000000000000000 -1.0000000000000000, 
0.0000000000000000 0.0000000000000000))'

            

      
      
      
    

        
True

So this looks like GEOS, not rgeos. I'll post to geos-devel for 
clarification.

Roger

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Fri, 25 Sep 2015, Vijay Lulla wrote:

        

            

      
      
      
    

        
Thanks for this, but are these the same - that is each geometry collection 
consists of two touching squares, but the second is offset =.1 North 
Eastwards? In any case, very grateful that you took a look!

Roger

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Edzer Pebesma
    

    

      
      #
    

  


  

      
Thanks Vijay; FYI, following OP's example, I see in PostGIS:

SELECT ST_AsText(ST_Intersection(
    'MULTIPOLYGON(
       ((0 0, 0 10, 10 10, 10 0, 0 0)),
       ((0 0, 10 0, 10 -10, 0 -10, 0 0)))'::geometry,
    'MULTIPOLYGON(
       ((0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5)),
       ((0.5 0.5, 10.5 0.5, 10.5 -10.5, 0.5 -10.5, 0.5 0.5)))'::geometry
    ));

giving:

ERROR:  Error performing intersection: TopologyException: Input geom 0
is invalid: Self-intersection at or near point 0 0 at 0 0

but

SELECT ST_AsText(ST_Intersection(
    'MULTIPOLYGON(
       ((0 0, 0 10, 10 10, 10 0, 0 0)),
       ((0 -1, 10 -1, 10 -10, 0 -10, 0 -1)))'::geometry,
    'MULTIPOLYGON(
       ((0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5)),
       ((0.5 0.0, 10.5 0.0, 10.5 -10.5, 0.5 -10.5, 0.5 0.0)))'::geometry
    ));

giving no error; when I try

SELECT ST_AsText(ST_Intersection(
    'POLYGON(
       (0 0, 0 10, 10 10, 10 0, 0 0),
       (0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry,
    'POLYGON(
       (0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5),
       (0.5 0.5, 10.5 0.5, 10.5 -10.5, 0.5 -10.5, 0.5 0.5))'::geometry
    ));

it gives the same error as the first example.

        
On 25/09/15 19:44, Roger Bivand wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Edzer Pebesma
    

    

      
      #
    

  


  

      
Trying to reduce the problem, it seems PostGIS thinks that polygons are
invalid for which R believes they are valid:

postgis=# select ST_IsValid('MULTIPOLYGON(
postgis'#        ((0 0, 0 10, 10 10, 10 0, 0 0)),
postgis'#        ((0 0, 10 0, 10 -10, 0 -10, 0 0)))'::geometry);
NOTICE:  Self-intersection at or near point 0 0
 st_isvalid
------------
 f
(1 row)

postgis=# select ST_IsValid('POLYGON(
postgis'#         (0 0, 0 10, 10 10, 10 0, 0 0),
postgis'#         (0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
NOTICE:  Self-intersection at or near point 0 0
 st_isvalid
------------
 f
(1 row)

but

            

      
      
      
    

        
[1] TRUE

Maybe gUnaryUnion-ing MultiPolygons before intersecting when byid =
FALSE is not such a bad idea after all?

        
On 25/09/15 21:55, Edzer Pebesma wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Vijay Lulla
    

    

      
      #
    

  


  

      
OOPS!  Sorry for my mistake.  I'm observing the same results that
Edzer has listed.  It has to do with validity
(http://www.postgis.net/docs/using_postgis_dbmanagement.html#OGC_Validity
) but I don't understand why the multipolygon listed in figure n
(similar to what we're examining here, I think) in section 4.3.5 of
postgis manual is invalid!

Please pardon the noise.

On Fri, Sep 25, 2015 at 3:55 PM, Edzer Pebesma

        
<edzer.pebesma at uni-muenster.de> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Edzer Pebesma
    

    

      
      #
    

  


  

      
On 25/09/15 23:30, Vijay Lulla wrote:

            

      
      
      
    

        
It says: "The boundaries of any two elements may touch, but only at a
finite number of POINTs." Points on the line that have both polygons in
common are ambiguous in terms of the DE-9IM relation with the polygon:
they are both inside the polygon as well as on the polygon boundary.

Still need to sort out why rgeos reports them as valid.

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Fri, 25 Sep 2015, Edzer Pebesma wrote:

        

            

      
      
      
    

        
Commited to R-forge as revision 505, seems to address the case, please 
check. The fix doesn't branch on geometry types, could those affected 
please check whether this leads to regressions?

Roger