[R-sig-Geo] gIntersection returns error "TopologyException: no outgoing dirEdge found at"
Jean-Luc Dupouey
dupouey at nancy.inra.fr
Fri Sep 25 13:16:58 CEST 2015
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,
>>
>>
>
More information about the R-sig-Geo
mailing list