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

[Roger Bivand]

On Fri, 25 Sep 2015, Roger Bivand wrote:

> On Fri, 25 Sep 2015, Jean-Luc Dupouey wrote:
>
...
>> 
>> 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))

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

> 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

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

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