gIntersection returns error "TopologyException: no outgoing dirEdge found at"

classic Classic list List threaded Threaded
15 messages Options
Reply | Threaded
Open this post in threaded view
|

gIntersection returns error "TopologyException: no outgoing dirEdge found at"

Jean-Luc DUPOUEY
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

It works with the option byid=TRUE.

But my question is why it does not work without? Is this behaviour
predictable?

I went through some previous related posts to the list, but could not
find the answer.

Thanks,

--
Jean-Luc Dupouey

INRA
Forest Ecology & Ecophysiology Unit
F-54280 Champenoux
France
mail: [hidden email]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: gIntersection returns error "TopologyException: no outgoing dirEdge found at"

Roger Bivand
Administrator
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,
>
>

--
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: [hidden email]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: gIntersection returns error "TopologyException: no outgoing dirEdge found at"

Jean-Luc DUPOUEY
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:

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)

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

Thanks in advance,

Jean-Luc Dupouey

INRA
Forest Ecology & Ecophysiology Unit
F-54280 Champenoux
France
mail: [hidden email]

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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: gIntersection returns error "TopologyException: no outgoing dirEdge found at"

Roger Bivand
Administrator
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: [hidden email]
>
> 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
> [hidden email]
> 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: [hidden email]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: gIntersection returns error "TopologyException: no outgoing dirEdge found at"

Roger Bivand
Administrator
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: [hidden email]
>>
>> 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
>> [hidden email]
>> 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: [hidden email]
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: gIntersection returns error "TopologyException: no outgoing dirEdge found at"

edzer


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: [hidden email]
>>>
>>> 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
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>>
>
>
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> 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


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

signature.asc (501 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: gIntersection returns error "TopologyException: no outgoing dirEdge found at"

Roger Bivand
Administrator
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: [hidden email]
>>>>
>>>> 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
>>>> [hidden email]
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>>>
>>
>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> 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: [hidden email]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: gIntersection returns error "TopologyException: no outgoing dirEdge found at"

Vijay Lulla
Seems to be working in PostGIS!  Below is from my PostGIS session:

gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 10, 10 10, 10 0, 0
0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
 st_intersects
---------------
 t
(1 row)


gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 10, 10 10, 10 0, 0
0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
                                  st_intersection
------------------------------------------------------------------------------------
 0102000000020000000000000000002440000000000000000000000000000000000000000000000000
(1 row)

gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 1, 1 1, 1 0, 0
0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
                                  st_intersection
------------------------------------------------------------------------------------
 010200000002000000000000000000F03F000000000000000000000000000000000000000000000000
(1 row)


gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 1, 1 1, 1 0, 0
0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
 st_intersects
---------------
 t
(1 row)

gisdb=# select postgis_version();
            postgis_version
---------------------------------------
 2.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1
(1 row)


gisdb=# select postgis_full_version();

  postgis_full_vers
ion
------------------------------------------------------------------------------------------
----------------------------------------------------------------------------
 POSTGIS="2.1.7 r13414" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel.
4.8.0, 6 March 2012" GDAL
="GDAL 1.11.1, released 2014/09/24" LIBXML="2.7.8" LIBJSON="UNKNOWN" RASTER
(1 row)

On Fri, Sep 25, 2015 at 1:02 PM, Roger Bivand <[hidden email]> 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).
>
> 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: [hidden email]
>>>>>
>>>>> 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
>>>>> [hidden email]
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>
>>>>
>>>>
>>>
>>>
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> [hidden email]
>>> 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: [hidden email]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: gIntersection returns error "TopologyException: no outgoing dirEdge found at"

Roger Bivand
Administrator
In reply to this post by Roger Bivand
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: [hidden email]
>>>>>
>>>>> 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
>>>>> [hidden email]
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>
>>>>
>>>>
>>>
>>>
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> [hidden email]
>>> 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: [hidden email]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: gIntersection returns error "TopologyException: no outgoing dirEdge found at"

Roger Bivand
Administrator
In reply to this post by Vijay Lulla
On Fri, 25 Sep 2015, Vijay Lulla wrote:

> Seems to be working in PostGIS!  Below is from my PostGIS session:
>
> gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 10, 10 10, 10 0, 0
> 0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
> st_intersects
> ---------------
> t
> (1 row)
>
>
> gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 10, 10 10, 10 0, 0
> 0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
>                                  st_intersection
> ------------------------------------------------------------------------------------
> 0102000000020000000000000000002440000000000000000000000000000000000000000000000000
> (1 row)
>
> gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 1, 1 1, 1 0, 0
> 0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
>                                  st_intersection
> ------------------------------------------------------------------------------------
> 010200000002000000000000000000F03F000000000000000000000000000000000000000000000000
> (1 row)
Thanks for this, but are these the same - that is each geometry collection
consists of two touching squares, but the second is offset =.1 North
Eastwards? In any case, very grateful that you took a look!

Roger

>
>
> gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 1, 1 1, 1 0, 0
> 0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
> st_intersects
> ---------------
> t
> (1 row)
>
> gisdb=# select postgis_version();
>            postgis_version
> ---------------------------------------
> 2.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1
> (1 row)
>
>
> gisdb=# select postgis_full_version();
>
>  postgis_full_vers
> ion
> ------------------------------------------------------------------------------------------
> ----------------------------------------------------------------------------
> POSTGIS="2.1.7 r13414" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel.
> 4.8.0, 6 March 2012" GDAL
> ="GDAL 1.11.1, released 2014/09/24" LIBXML="2.7.8" LIBJSON="UNKNOWN" RASTER
> (1 row)
>
> On Fri, Sep 25, 2015 at 1:02 PM, Roger Bivand <[hidden email]> 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).
>>
>> 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: [hidden email]
>>>>>>
>>>>>> 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
>>>>>> [hidden email]
>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> [hidden email]
>>>> 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: [hidden email]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> 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: [hidden email]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: gIntersection returns error "TopologyException: no outgoing dirEdge found at"

edzer
Thanks Vijay; FYI, following OP's example, I see in PostGIS:

SELECT ST_AsText(ST_Intersection(
    'MULTIPOLYGON(
       ((0 0, 0 10, 10 10, 10 0, 0 0)),
       ((0 0, 10 0, 10 -10, 0 -10, 0 0)))'::geometry,
    'MULTIPOLYGON(
       ((0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5)),
       ((0.5 0.5, 10.5 0.5, 10.5 -10.5, 0.5 -10.5, 0.5 0.5)))'::geometry
    ));

giving:

ERROR:  Error performing intersection: TopologyException: Input geom 0
is invalid: Self-intersection at or near point 0 0 at 0 0

but

SELECT ST_AsText(ST_Intersection(
    'MULTIPOLYGON(
       ((0 0, 0 10, 10 10, 10 0, 0 0)),
       ((0 -1, 10 -1, 10 -10, 0 -10, 0 -1)))'::geometry,
    'MULTIPOLYGON(
       ((0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5)),
       ((0.5 0.0, 10.5 0.0, 10.5 -10.5, 0.5 -10.5, 0.5 0.0)))'::geometry
    ));

giving no error; when I try

SELECT ST_AsText(ST_Intersection(
    'POLYGON(
       (0 0, 0 10, 10 10, 10 0, 0 0),
       (0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry,
    'POLYGON(
       (0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5),
       (0.5 0.5, 10.5 0.5, 10.5 -10.5, 0.5 -10.5, 0.5 0.5))'::geometry
    ));

it gives the same error as the first example.



On 25/09/15 19:44, Roger Bivand wrote:

> On Fri, 25 Sep 2015, Vijay Lulla wrote:
>
>> Seems to be working in PostGIS!  Below is from my PostGIS session:
>>
>> gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 10, 10 10, 10 0, 0
>> 0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
>> st_intersects
>> ---------------
>> t
>> (1 row)
>>
>>
>> gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 10, 10 10, 10 0, 0
>> 0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
>>                                  st_intersection
>> ------------------------------------------------------------------------------------
>>
>> 0102000000020000000000000000002440000000000000000000000000000000000000000000000000
>>
>> (1 row)
>>
>> gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 1, 1 1, 1 0, 0
>> 0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
>>                                  st_intersection
>> ------------------------------------------------------------------------------------
>>
>> 010200000002000000000000000000F03F000000000000000000000000000000000000000000000000
>>
>> (1 row)
>
> Thanks for this, but are these the same - that is each geometry
> collection consists of two touching squares, but the second is offset
> =.1 North Eastwards? In any case, very grateful that you took a look!
>
> Roger
>
>>
>>
>> gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 1, 1 1, 1 0, 0
>> 0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
>> st_intersects
>> ---------------
>> t
>> (1 row)
>>
>> gisdb=# select postgis_version();
>>            postgis_version
>> ---------------------------------------
>> 2.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1
>> (1 row)
>>
>>
>> gisdb=# select postgis_full_version();
>>
>>  postgis_full_vers
>> ion
>> ------------------------------------------------------------------------------------------
>>
>> ----------------------------------------------------------------------------
>>
>> POSTGIS="2.1.7 r13414" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel.
>> 4.8.0, 6 March 2012" GDAL
>> ="GDAL 1.11.1, released 2014/09/24" LIBXML="2.7.8" LIBJSON="UNKNOWN"
>> RASTER
>> (1 row)
>>
>> On Fri, Sep 25, 2015 at 1:02 PM, Roger Bivand <[hidden email]>
>> 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).
>>>
>>> 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: [hidden email]
>>>>>>>
>>>>>>> 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
>>>>>>> [hidden email]
>>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>>
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-Geo mailing list
>>>>> [hidden email]
>>>>> 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: [hidden email]
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> [hidden email]
>>> 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


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

signature.asc (501 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: gIntersection returns error "TopologyException: no outgoing dirEdge found at"

edzer
Trying to reduce the problem, it seems PostGIS thinks that polygons are
invalid for which R believes they are valid:

postgis=# select ST_IsValid('MULTIPOLYGON(
postgis'#        ((0 0, 0 10, 10 10, 10 0, 0 0)),
postgis'#        ((0 0, 10 0, 10 -10, 0 -10, 0 0)))'::geometry);
NOTICE:  Self-intersection at or near point 0 0
 st_isvalid
------------
 f
(1 row)

postgis=# select ST_IsValid('POLYGON(
postgis'#         (0 0, 0 10, 10 10, 10 0, 0 0),
postgis'#         (0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
NOTICE:  Self-intersection at or near point 0 0
 st_isvalid
------------
 f
(1 row)

but

> 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))
>
> Pols1=Polygons(list(Polygon(Pol1)),"Pols1")
> Pols2=Polygons(list(Polygon(Pol2)),"Pols2")
> MyLay=SpatialPolygons(list(Pols1,Pols2))
> gIsValid(MyLay)
[1] TRUE

Maybe gUnaryUnion-ing MultiPolygons before intersecting when byid =
FALSE is not such a bad idea after all?


On 25/09/15 21:55, Edzer Pebesma wrote:

> Thanks Vijay; FYI, following OP's example, I see in PostGIS:
>
> SELECT ST_AsText(ST_Intersection(
>     'MULTIPOLYGON(
>        ((0 0, 0 10, 10 10, 10 0, 0 0)),
>        ((0 0, 10 0, 10 -10, 0 -10, 0 0)))'::geometry,
>     'MULTIPOLYGON(
>        ((0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5)),
>        ((0.5 0.5, 10.5 0.5, 10.5 -10.5, 0.5 -10.5, 0.5 0.5)))'::geometry
>     ));
>
> giving:
>
> ERROR:  Error performing intersection: TopologyException: Input geom 0
> is invalid: Self-intersection at or near point 0 0 at 0 0
>
> but
>
> SELECT ST_AsText(ST_Intersection(
>     'MULTIPOLYGON(
>        ((0 0, 0 10, 10 10, 10 0, 0 0)),
>        ((0 -1, 10 -1, 10 -10, 0 -10, 0 -1)))'::geometry,
>     'MULTIPOLYGON(
>        ((0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5)),
>        ((0.5 0.0, 10.5 0.0, 10.5 -10.5, 0.5 -10.5, 0.5 0.0)))'::geometry
>     ));
>
> giving no error; when I try
>
> SELECT ST_AsText(ST_Intersection(
>     'POLYGON(
>        (0 0, 0 10, 10 10, 10 0, 0 0),
>        (0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry,
>     'POLYGON(
>        (0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5),
>        (0.5 0.5, 10.5 0.5, 10.5 -10.5, 0.5 -10.5, 0.5 0.5))'::geometry
>     ));
>
> it gives the same error as the first example.
>
>
>
> On 25/09/15 19:44, Roger Bivand wrote:
>> On Fri, 25 Sep 2015, Vijay Lulla wrote:
>>
>>> Seems to be working in PostGIS!  Below is from my PostGIS session:
>>>
>>> gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 10, 10 10, 10 0, 0
>>> 0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
>>> st_intersects
>>> ---------------
>>> t
>>> (1 row)
>>>
>>>
>>> gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 10, 10 10, 10 0, 0
>>> 0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
>>>                                  st_intersection
>>> ------------------------------------------------------------------------------------
>>>
>>> 0102000000020000000000000000002440000000000000000000000000000000000000000000000000
>>>
>>> (1 row)
>>>
>>> gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 1, 1 1, 1 0, 0
>>> 0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
>>>                                  st_intersection
>>> ------------------------------------------------------------------------------------
>>>
>>> 010200000002000000000000000000F03F000000000000000000000000000000000000000000000000
>>>
>>> (1 row)
>>
>> Thanks for this, but are these the same - that is each geometry
>> collection consists of two touching squares, but the second is offset
>> =.1 North Eastwards? In any case, very grateful that you took a look!
>>
>> Roger
>>
>>>
>>>
>>> gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 1, 1 1, 1 0, 0
>>> 0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
>>> st_intersects
>>> ---------------
>>> t
>>> (1 row)
>>>
>>> gisdb=# select postgis_version();
>>>            postgis_version
>>> ---------------------------------------
>>> 2.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1
>>> (1 row)
>>>
>>>
>>> gisdb=# select postgis_full_version();
>>>
>>>  postgis_full_vers
>>> ion
>>> ------------------------------------------------------------------------------------------
>>>
>>> ----------------------------------------------------------------------------
>>>
>>> POSTGIS="2.1.7 r13414" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel.
>>> 4.8.0, 6 March 2012" GDAL
>>> ="GDAL 1.11.1, released 2014/09/24" LIBXML="2.7.8" LIBJSON="UNKNOWN"
>>> RASTER
>>> (1 row)
>>>
>>> On Fri, Sep 25, 2015 at 1:02 PM, Roger Bivand <[hidden email]>
>>> 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).
>>>>
>>>> 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: [hidden email]
>>>>>>>>
>>>>>>>> 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
>>>>>>>> [hidden email]
>>>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> _______________________________________________
>>>>>> R-sig-Geo mailing list
>>>>>> [hidden email]
>>>>>> 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: [hidden email]
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> [hidden email]
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>>
>
>
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> 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


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

signature.asc (501 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: gIntersection returns error "TopologyException: no outgoing dirEdge found at"

Vijay Lulla
In reply to this post by edzer
OOPS!  Sorry for my mistake.  I'm observing the same results that
Edzer has listed.  It has to do with validity
(http://www.postgis.net/docs/using_postgis_dbmanagement.html#OGC_Validity
) but I don't understand why the multipolygon listed in figure n
(similar to what we're examining here, I think) in section 4.3.5 of
postgis manual is invalid!

Please pardon the noise.

On Fri, Sep 25, 2015 at 3:55 PM, Edzer Pebesma
<[hidden email]> wrote:

> Thanks Vijay; FYI, following OP's example, I see in PostGIS:
>
> SELECT ST_AsText(ST_Intersection(
>     'MULTIPOLYGON(
>        ((0 0, 0 10, 10 10, 10 0, 0 0)),
>        ((0 0, 10 0, 10 -10, 0 -10, 0 0)))'::geometry,
>     'MULTIPOLYGON(
>        ((0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5)),
>        ((0.5 0.5, 10.5 0.5, 10.5 -10.5, 0.5 -10.5, 0.5 0.5)))'::geometry
>     ));
>
> giving:
>
> ERROR:  Error performing intersection: TopologyException: Input geom 0
> is invalid: Self-intersection at or near point 0 0 at 0 0
>
> but
>
> SELECT ST_AsText(ST_Intersection(
>     'MULTIPOLYGON(
>        ((0 0, 0 10, 10 10, 10 0, 0 0)),
>        ((0 -1, 10 -1, 10 -10, 0 -10, 0 -1)))'::geometry,
>     'MULTIPOLYGON(
>        ((0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5)),
>        ((0.5 0.0, 10.5 0.0, 10.5 -10.5, 0.5 -10.5, 0.5 0.0)))'::geometry
>     ));
>
> giving no error; when I try
>
> SELECT ST_AsText(ST_Intersection(
>     'POLYGON(
>        (0 0, 0 10, 10 10, 10 0, 0 0),
>        (0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry,
>     'POLYGON(
>        (0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5),
>        (0.5 0.5, 10.5 0.5, 10.5 -10.5, 0.5 -10.5, 0.5 0.5))'::geometry
>     ));
>
> it gives the same error as the first example.
>
>
>
> On 25/09/15 19:44, Roger Bivand wrote:
>> On Fri, 25 Sep 2015, Vijay Lulla wrote:
>>
>>> Seems to be working in PostGIS!  Below is from my PostGIS session:
>>>
>>> gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 10, 10 10, 10 0, 0
>>> 0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
>>> st_intersects
>>> ---------------
>>> t
>>> (1 row)
>>>
>>>
>>> gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 10, 10 10, 10 0, 0
>>> 0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
>>>                                  st_intersection
>>> ------------------------------------------------------------------------------------
>>>
>>> 0102000000020000000000000000002440000000000000000000000000000000000000000000000000
>>>
>>> (1 row)
>>>
>>> gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 1, 1 1, 1 0, 0
>>> 0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
>>>                                  st_intersection
>>> ------------------------------------------------------------------------------------
>>>
>>> 010200000002000000000000000000F03F000000000000000000000000000000000000000000000000
>>>
>>> (1 row)
>>
>> Thanks for this, but are these the same - that is each geometry
>> collection consists of two touching squares, but the second is offset
>> =.1 North Eastwards? In any case, very grateful that you took a look!
>>
>> Roger
>>
>>>
>>>
>>> gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 1, 1 1, 1 0, 0
>>> 0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
>>> st_intersects
>>> ---------------
>>> t
>>> (1 row)
>>>
>>> gisdb=# select postgis_version();
>>>            postgis_version
>>> ---------------------------------------
>>> 2.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1
>>> (1 row)
>>>
>>>
>>> gisdb=# select postgis_full_version();
>>>
>>>  postgis_full_vers
>>> ion
>>> ------------------------------------------------------------------------------------------
>>>
>>> ----------------------------------------------------------------------------
>>>
>>> POSTGIS="2.1.7 r13414" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel.
>>> 4.8.0, 6 March 2012" GDAL
>>> ="GDAL 1.11.1, released 2014/09/24" LIBXML="2.7.8" LIBJSON="UNKNOWN"
>>> RASTER
>>> (1 row)
>>>
>>> On Fri, Sep 25, 2015 at 1:02 PM, Roger Bivand <[hidden email]>
>>> 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).
>>>>
>>>> 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: [hidden email]
>>>>>>>>
>>>>>>>> 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
>>>>>>>> [hidden email]
>>>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> _______________________________________________
>>>>>> R-sig-Geo mailing list
>>>>>> [hidden email]
>>>>>> 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: [hidden email]
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> [hidden email]
>>>> 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
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: gIntersection returns error "TopologyException: no outgoing dirEdge found at"

edzer


On 25/09/15 23:30, Vijay Lulla wrote:
> OOPS!  Sorry for my mistake.  I'm observing the same results that
> Edzer has listed.  It has to do with validity
> (http://www.postgis.net/docs/using_postgis_dbmanagement.html#OGC_Validity
> ) but I don't understand why the multipolygon listed in figure n
> (similar to what we're examining here, I think) in section 4.3.5 of
> postgis manual is invalid!

It says: "The boundaries of any two elements may touch, but only at a
finite number of POINTs." Points on the line that have both polygons in
common are ambiguous in terms of the DE-9IM relation with the polygon:
they are both inside the polygon as well as on the polygon boundary.

Still need to sort out why rgeos reports them as valid.

>
> Please pardon the noise.
>
> On Fri, Sep 25, 2015 at 3:55 PM, Edzer Pebesma
> <[hidden email]> wrote:
>> Thanks Vijay; FYI, following OP's example, I see in PostGIS:
>>
>> SELECT ST_AsText(ST_Intersection(
>>     'MULTIPOLYGON(
>>        ((0 0, 0 10, 10 10, 10 0, 0 0)),
>>        ((0 0, 10 0, 10 -10, 0 -10, 0 0)))'::geometry,
>>     'MULTIPOLYGON(
>>        ((0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5)),
>>        ((0.5 0.5, 10.5 0.5, 10.5 -10.5, 0.5 -10.5, 0.5 0.5)))'::geometry
>>     ));
>>
>> giving:
>>
>> ERROR:  Error performing intersection: TopologyException: Input geom 0
>> is invalid: Self-intersection at or near point 0 0 at 0 0
>>
>> but
>>
>> SELECT ST_AsText(ST_Intersection(
>>     'MULTIPOLYGON(
>>        ((0 0, 0 10, 10 10, 10 0, 0 0)),
>>        ((0 -1, 10 -1, 10 -10, 0 -10, 0 -1)))'::geometry,
>>     'MULTIPOLYGON(
>>        ((0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5)),
>>        ((0.5 0.0, 10.5 0.0, 10.5 -10.5, 0.5 -10.5, 0.5 0.0)))'::geometry
>>     ));
>>
>> giving no error; when I try
>>
>> SELECT ST_AsText(ST_Intersection(
>>     'POLYGON(
>>        (0 0, 0 10, 10 10, 10 0, 0 0),
>>        (0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry,
>>     'POLYGON(
>>        (0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5),
>>        (0.5 0.5, 10.5 0.5, 10.5 -10.5, 0.5 -10.5, 0.5 0.5))'::geometry
>>     ));
>>
>> it gives the same error as the first example.
>>
>>
>>
>> On 25/09/15 19:44, Roger Bivand wrote:
>>> On Fri, 25 Sep 2015, Vijay Lulla wrote:
>>>
>>>> Seems to be working in PostGIS!  Below is from my PostGIS session:
>>>>
>>>> gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 10, 10 10, 10 0, 0
>>>> 0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
>>>> st_intersects
>>>> ---------------
>>>> t
>>>> (1 row)
>>>>
>>>>
>>>> gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 10, 10 10, 10 0, 0
>>>> 0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
>>>>                                  st_intersection
>>>> ------------------------------------------------------------------------------------
>>>>
>>>> 0102000000020000000000000000002440000000000000000000000000000000000000000000000000
>>>>
>>>> (1 row)
>>>>
>>>> gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 1, 1 1, 1 0, 0
>>>> 0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
>>>>                                  st_intersection
>>>> ------------------------------------------------------------------------------------
>>>>
>>>> 010200000002000000000000000000F03F000000000000000000000000000000000000000000000000
>>>>
>>>> (1 row)
>>>
>>> Thanks for this, but are these the same - that is each geometry
>>> collection consists of two touching squares, but the second is offset
>>> =.1 North Eastwards? In any case, very grateful that you took a look!
>>>
>>> Roger
>>>
>>>>
>>>>
>>>> gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 1, 1 1, 1 0, 0
>>>> 0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
>>>> st_intersects
>>>> ---------------
>>>> t
>>>> (1 row)
>>>>
>>>> gisdb=# select postgis_version();
>>>>            postgis_version
>>>> ---------------------------------------
>>>> 2.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1
>>>> (1 row)
>>>>
>>>>
>>>> gisdb=# select postgis_full_version();
>>>>
>>>>  postgis_full_vers
>>>> ion
>>>> ------------------------------------------------------------------------------------------
>>>>
>>>> ----------------------------------------------------------------------------
>>>>
>>>> POSTGIS="2.1.7 r13414" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel.
>>>> 4.8.0, 6 March 2012" GDAL
>>>> ="GDAL 1.11.1, released 2014/09/24" LIBXML="2.7.8" LIBJSON="UNKNOWN"
>>>> RASTER
>>>> (1 row)
>>>>
>>>> On Fri, Sep 25, 2015 at 1:02 PM, Roger Bivand <[hidden email]>
>>>> 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).
>>>>>
>>>>> 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: [hidden email]
>>>>>>>>>
>>>>>>>>> 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
>>>>>>>>> [hidden email]
>>>>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> R-sig-Geo mailing list
>>>>>>> [hidden email]
>>>>>>> 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: [hidden email]
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-Geo mailing list
>>>>> [hidden email]
>>>>> 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
>>
--
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


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

signature.asc (501 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: gIntersection returns error "TopologyException: no outgoing dirEdge found at"

Roger Bivand
Administrator
In reply to this post by edzer
On Fri, 25 Sep 2015, Edzer Pebesma wrote:

> Trying to reduce the problem, it seems PostGIS thinks that polygons are
> invalid for which R believes they are valid:
>
> postgis=# select ST_IsValid('MULTIPOLYGON(
> postgis'#        ((0 0, 0 10, 10 10, 10 0, 0 0)),
> postgis'#        ((0 0, 10 0, 10 -10, 0 -10, 0 0)))'::geometry);
> NOTICE:  Self-intersection at or near point 0 0
> st_isvalid
> ------------
> f
> (1 row)
>
> postgis=# select ST_IsValid('POLYGON(
> postgis'#         (0 0, 0 10, 10 10, 10 0, 0 0),
> postgis'#         (0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
> NOTICE:  Self-intersection at or near point 0 0
> st_isvalid
> ------------
> f
> (1 row)
>
> but
>
>> 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))
>>
>> Pols1=Polygons(list(Polygon(Pol1)),"Pols1")
>> Pols2=Polygons(list(Polygon(Pol2)),"Pols2")
>> MyLay=SpatialPolygons(list(Pols1,Pols2))
>> gIsValid(MyLay)
> [1] TRUE
>
> Maybe gUnaryUnion-ing MultiPolygons before intersecting when byid =
> FALSE is not such a bad idea after all?
>
Commited to R-forge as revision 505, seems to address the case, please
check. The fix doesn't branch on geometry types, could those affected
please check whether this leads to regressions?

Roger

>
> On 25/09/15 21:55, Edzer Pebesma wrote:
>> Thanks Vijay; FYI, following OP's example, I see in PostGIS:
>>
>> SELECT ST_AsText(ST_Intersection(
>>     'MULTIPOLYGON(
>>        ((0 0, 0 10, 10 10, 10 0, 0 0)),
>>        ((0 0, 10 0, 10 -10, 0 -10, 0 0)))'::geometry,
>>     'MULTIPOLYGON(
>>        ((0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5)),
>>        ((0.5 0.5, 10.5 0.5, 10.5 -10.5, 0.5 -10.5, 0.5 0.5)))'::geometry
>>     ));
>>
>> giving:
>>
>> ERROR:  Error performing intersection: TopologyException: Input geom 0
>> is invalid: Self-intersection at or near point 0 0 at 0 0
>>
>> but
>>
>> SELECT ST_AsText(ST_Intersection(
>>     'MULTIPOLYGON(
>>        ((0 0, 0 10, 10 10, 10 0, 0 0)),
>>        ((0 -1, 10 -1, 10 -10, 0 -10, 0 -1)))'::geometry,
>>     'MULTIPOLYGON(
>>        ((0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5)),
>>        ((0.5 0.0, 10.5 0.0, 10.5 -10.5, 0.5 -10.5, 0.5 0.0)))'::geometry
>>     ));
>>
>> giving no error; when I try
>>
>> SELECT ST_AsText(ST_Intersection(
>>     'POLYGON(
>>        (0 0, 0 10, 10 10, 10 0, 0 0),
>>        (0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry,
>>     'POLYGON(
>>        (0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5),
>>        (0.5 0.5, 10.5 0.5, 10.5 -10.5, 0.5 -10.5, 0.5 0.5))'::geometry
>>     ));
>>
>> it gives the same error as the first example.
>>
>>
>>
>> On 25/09/15 19:44, Roger Bivand wrote:
>>> On Fri, 25 Sep 2015, Vijay Lulla wrote:
>>>
>>>> Seems to be working in PostGIS!  Below is from my PostGIS session:
>>>>
>>>> gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 10, 10 10, 10 0, 0
>>>> 0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
>>>> st_intersects
>>>> ---------------
>>>> t
>>>> (1 row)
>>>>
>>>>
>>>> gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 10, 10 10, 10 0, 0
>>>> 0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
>>>>                                  st_intersection
>>>> ------------------------------------------------------------------------------------
>>>>
>>>> 0102000000020000000000000000002440000000000000000000000000000000000000000000000000
>>>>
>>>> (1 row)
>>>>
>>>> gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 1, 1 1, 1 0, 0
>>>> 0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
>>>>                                  st_intersection
>>>> ------------------------------------------------------------------------------------
>>>>
>>>> 010200000002000000000000000000F03F000000000000000000000000000000000000000000000000
>>>>
>>>> (1 row)
>>>
>>> Thanks for this, but are these the same - that is each geometry
>>> collection consists of two touching squares, but the second is offset
>>> =.1 North Eastwards? In any case, very grateful that you took a look!
>>>
>>> Roger
>>>
>>>>
>>>>
>>>> gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 1, 1 1, 1 0, 0
>>>> 0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
>>>> st_intersects
>>>> ---------------
>>>> t
>>>> (1 row)
>>>>
>>>> gisdb=# select postgis_version();
>>>>            postgis_version
>>>> ---------------------------------------
>>>> 2.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1
>>>> (1 row)
>>>>
>>>>
>>>> gisdb=# select postgis_full_version();
>>>>
>>>>  postgis_full_vers
>>>> ion
>>>> ------------------------------------------------------------------------------------------
>>>>
>>>> ----------------------------------------------------------------------------
>>>>
>>>> POSTGIS="2.1.7 r13414" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel.
>>>> 4.8.0, 6 March 2012" GDAL
>>>> ="GDAL 1.11.1, released 2014/09/24" LIBXML="2.7.8" LIBJSON="UNKNOWN"
>>>> RASTER
>>>> (1 row)
>>>>
>>>> On Fri, Sep 25, 2015 at 1:02 PM, Roger Bivand <[hidden email]>
>>>> 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).
>>>>>
>>>>> 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: [hidden email]
>>>>>>>>>
>>>>>>>>> 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
>>>>>>>>> [hidden email]
>>>>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> R-sig-Geo mailing list
>>>>>>> [hidden email]
>>>>>>> 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: [hidden email]
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-Geo mailing list
>>>>> [hidden email]
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>
>>>>
>>>
>>
>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> 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: [hidden email]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway