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

[Edzer Pebesma]

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?

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

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi),  University of Münster,
Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
Journal of Statistical Software:   http://www.jstatsoft.org/
Computers & Geosciences:   http://elsevier.com/locate/cageo/
Spatial Statistics Society http://www.spatialstatistics.info

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 490 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20150925/64eb6168/attachment.bin>