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

Roger Bivand Roger.Bivand at nhh.no
Fri Sep 25 19:29:41 CEST 2015 

On Fri, 25 Sep 2015, Roger Bivand wrote:

> On Fri, 25 Sep 2015, Edzer Pebesma wrote:
>
>> 
>> 
>> On 25/09/15 14:36, Roger Bivand wrote:
>>> 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?
>> 
>> Has someone tried to reproduce this with GEOS, without R, e.g. by a
>> PostGIS query?
>
> Do you have a running PostGIS instance - I don't. I did try v.overlay in 
> GRASS, but it isn't GEOS-based, and returns the correct answer for 
> byid=c(TRUE, TRUE) so like rgeos in the same setting. I'm re-installing GEOS 
> with Python, but someone else should check (my GEOS 3.5.0).

With Shapely, I see (values output from rgeos::writeWKT):

from shapely.wkt import dumps, loads
>>> MyLay = loads("GEOMETRYCOLLECTION (POLYGON ((0.0000000000000000 
0.0000000000000000, 0.0000000000000000 1.0000000000000000, 
1.0000000000000000 1.0000000000000000, 1.0000000000000000 
0.0000000000000000, 0.0000000000000000 0.0000000000000000)), POLYGON 
((0.0000000000000000 0.0000000000000000, 1.0000000000000000 
0.0000000000000000, 1.0000000000000000 -1.0000000000000000, 
0.0000000000000000 -1.0000000000000000, 0.0000000000000000 
0.0000000000000000)))")
>>> MyLayl = loads("GEOMETRYCOLLECTION (POLYGON ((0.1000000000000000 
0.1000000000000000, 0.1000000000000000 1.1000000000000001, 
1.1000000000000001 1.1000000000000001, 1.1000000000000001 
0.1000000000000000, 0.1000000000000000 0.1000000000000000)), POLYGON 
((0.1000000000000000 0.1000000000000000, 1.1000000000000001 
0.1000000000000000, 1.1000000000000001 -0.9000000000000000, 
0.1000000000000000 -0.9000000000000000, 0.1000000000000000 
0.1000000000000000)))")
>>> merged = MyLay.intersection(MyLayl)
>>> dumps(merged)
'POLYGON ((0.0000000000000000 0.0000000000000000, 0.0000000000000000 
1.0000000000000000, 1.0000000000000000 1.0000000000000000, 
1.0000000000000000 0.0000000000000000, 1.0000000000000000 
-1.0000000000000000, 0.0000000000000000 -1.0000000000000000, 
0.0000000000000000 0.0000000000000000))'
>>> merged.equals(MyLay)
True

So this looks like GEOS, not rgeos. I'll post to geos-devel for 
clarification.

Roger

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

More information about the R-sig-Geo mailing list