[R-sig-Geo] gIntersection returns error "TopologyException: no outgoing dirEdge found at"

[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 :
> 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,
>>
>>
>