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

Fri Sep 25 19:02:15 CEST 2015

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

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