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

[Roger Bivand]

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

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

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.

>
> 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)

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.

>
> #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.

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

> 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?
>

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

> 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,
>>> 
>>> 
>> 
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: Roger.Bivand at nhh.no