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

Fri, Sep 25, 2015 4:16 AM

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 :

On Wed, 23 Sep 2015, Jean-Luc Dupouey wrote:

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

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.

It works with the option byid=TRUE.

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

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

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

Thanks,

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

Thread (15 messages)