

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,

JeanLuc Dupouey
INRA
Forest Ecology & Ecophysiology Unit
F54280 Champenoux
France
mail: [hidden email]
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo

Administrator

On Wed, 23 Sep 2015, JeanLuc 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, N5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
email: [hidden email]
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N5045 Bergen, Norway


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.312 and sp_1.20.
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,
JeanLuc Dupouey
INRA
Forest Ecology & Ecophysiology Unit
F54280 Champenoux
France
mail: [hidden email]
Le 23/09/2015 19:39, Roger Bivand a écrit :
> On Wed, 23 Sep 2015, JeanLuc 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,
>>
>>
>
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo

Administrator

On Fri, 25 Sep 2015, JeanLuc 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.312 and sp_1.20.
>
> 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,
>
> JeanLuc Dupouey
>
> INRA
> Forest Ecology & Ecophysiology Unit
> F54280 Champenoux
> France
> mail: [hidden email]
>
> Le 23/09/2015 19:39, Roger Bivand a écrit :
>> On Wed, 23 Sep 2015, JeanLuc 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,
>>>
>>>
>>
>
> _______________________________________________
> RsigGeo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/rsiggeo>

Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
email: [hidden email]
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N5045 Bergen, Norway

Administrator

On Fri, 25 Sep 2015, Roger Bivand wrote:
> On Fri, 25 Sep 2015, JeanLuc 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.312 and sp_1.20.
>>
>> 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,
>>
>> JeanLuc Dupouey
>>
>> INRA
>> Forest Ecology & Ecophysiology Unit
>> F54280 Champenoux
>> France
>> mail: [hidden email]
>>
>> Le 23/09/2015 19:39, Roger Bivand a écrit :
>>> On Wed, 23 Sep 2015, JeanLuc 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,
>>>>
>>>>
>>>
>>
>> _______________________________________________
>> RsigGeo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>
>
>

Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
email: [hidden email]_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N5045 Bergen, Norway


On 25/09/15 14:36, Roger Bivand wrote:
> On Fri, 25 Sep 2015, Roger Bivand wrote:
>
>> On Fri, 25 Sep 2015, JeanLuc 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.312 and sp_1.20.
>>>
>>> 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,
>>>
>>> JeanLuc Dupouey
>>>
>>> INRA
>>> Forest Ecology & Ecophysiology Unit
>>> F54280 Champenoux
>>> France
>>> mail: [hidden email]
>>>
>>> Le 23/09/2015 19:39, Roger Bivand a écrit :
>>>> On Wed, 23 Sep 2015, JeanLuc 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,
>>>>>
>>>>>
>>>>
>>>
>>> _______________________________________________
>>> RsigGeo mailing list
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>
>>
>>
>
>
>
> _______________________________________________
> RsigGeo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/rsiggeo>

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_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo

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, JeanLuc 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 GEOSbased, and returns the correct answer for
byid=c(TRUE, TRUE) so like rgeos in the same setting. I'm reinstalling
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.312 and sp_1.20.
>>>>
>>>> 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,
>>>>
>>>> JeanLuc Dupouey
>>>>
>>>> INRA
>>>> Forest Ecology & Ecophysiology Unit
>>>> F54280 Champenoux
>>>> France
>>>> mail: [hidden email]
>>>>
>>>> Le 23/09/2015 19:39, Roger Bivand a écrit :
>>>>> On Wed, 23 Sep 2015, JeanLuc 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,
>>>>>>
>>>>>>
>>>>>
>>>>
>>>> _______________________________________________
>>>> RsigGeo mailing list
>>>> [hidden email]
>>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>>
>>>
>>>
>>
>>
>>
>> _______________________________________________
>> RsigGeo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>
>
>

Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
email: [hidden email]
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N5045 Bergen, Norway


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.2CAPI1.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, JeanLuc 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 GEOSbased, and returns the correct answer for
> byid=c(TRUE, TRUE) so like rgeos in the same setting. I'm reinstalling 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.312 and sp_1.20.
>>>>>
>>>>> 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,
>>>>>
>>>>> JeanLuc Dupouey
>>>>>
>>>>> INRA
>>>>> Forest Ecology & Ecophysiology Unit
>>>>> F54280 Champenoux
>>>>> France
>>>>> mail: [hidden email]
>>>>>
>>>>> Le 23/09/2015 19:39, Roger Bivand a écrit :
>>>>>>
>>>>>> On Wed, 23 Sep 2015, JeanLuc 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,
>>>>>>>
>>>>>>>
>>>>>>
>>>>>
>>>>> _______________________________________________
>>>>> RsigGeo mailing list
>>>>> [hidden email]
>>>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>>>
>>>>
>>>>
>>>
>>>
>>>
>>> _______________________________________________
>>> RsigGeo mailing list
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>
>>
>>
>
> 
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N5045 Bergen, Norway.
> voice: +47 55 95 93 55; fax +47 55 95 91 00
> email: [hidden email]
>
> _______________________________________________
> RsigGeo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/rsiggeo>
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo

Administrator

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, JeanLuc 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 GEOSbased, and returns the correct answer for
> byid=c(TRUE, TRUE) so like rgeos in the same setting. I'm reinstalling 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 geosdevel 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.312 and sp_1.20.
>>>>>
>>>>> 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,
>>>>>
>>>>> JeanLuc Dupouey
>>>>>
>>>>> INRA
>>>>> Forest Ecology & Ecophysiology Unit
>>>>> F54280 Champenoux
>>>>> France
>>>>> mail: [hidden email]
>>>>>
>>>>> Le 23/09/2015 19:39, Roger Bivand a écrit :
>>>>>> On Wed, 23 Sep 2015, JeanLuc 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,
>>>>>>>
>>>>>>>
>>>>>>
>>>>>
>>>>> _______________________________________________
>>>>> RsigGeo mailing list
>>>>> [hidden email]
>>>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>>>
>>>>
>>>>
>>>
>>>
>>>
>>> _______________________________________________
>>> RsigGeo mailing list
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>
>>
>>
>
>

Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
email: [hidden email]
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N5045 Bergen, Norway

Administrator

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.2CAPI1.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, JeanLuc 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 GEOSbased, and returns the correct answer for
>> byid=c(TRUE, TRUE) so like rgeos in the same setting. I'm reinstalling 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.312 and sp_1.20.
>>>>>>
>>>>>> 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,
>>>>>>
>>>>>> JeanLuc Dupouey
>>>>>>
>>>>>> INRA
>>>>>> Forest Ecology & Ecophysiology Unit
>>>>>> F54280 Champenoux
>>>>>> France
>>>>>> mail: [hidden email]
>>>>>>
>>>>>> Le 23/09/2015 19:39, Roger Bivand a écrit :
>>>>>>>
>>>>>>> On Wed, 23 Sep 2015, JeanLuc 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,
>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>
>>>>>> _______________________________________________
>>>>>> RsigGeo mailing list
>>>>>> [hidden email]
>>>>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>>
>>>> _______________________________________________
>>>> RsigGeo mailing list
>>>> [hidden email]
>>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>>
>>>
>>>
>>
>> 
>> Roger Bivand
>> Department of Economics, Norwegian School of Economics,
>> Helleveien 30, N5045 Bergen, Norway.
>> voice: +47 55 95 93 55; fax +47 55 95 91 00
>> email: [hidden email]
>>
>> _______________________________________________
>> RsigGeo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>
>

Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
email: [hidden email]
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N5045 Bergen, Norway


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: Selfintersection 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.2CAPI1.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, JeanLuc 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 GEOSbased, and returns the correct answer for
>>> byid=c(TRUE, TRUE) so like rgeos in the same setting. I'm
>>> reinstalling 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.312 and sp_1.20.
>>>>>>>
>>>>>>> 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,
>>>>>>>
>>>>>>> JeanLuc Dupouey
>>>>>>>
>>>>>>> INRA
>>>>>>> Forest Ecology & Ecophysiology Unit
>>>>>>> F54280 Champenoux
>>>>>>> France
>>>>>>> mail: [hidden email]
>>>>>>>
>>>>>>> Le 23/09/2015 19:39, Roger Bivand a écrit :
>>>>>>>>
>>>>>>>> On Wed, 23 Sep 2015, JeanLuc 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,
>>>>>>>>>
>>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> RsigGeo mailing list
>>>>>>> [hidden email]
>>>>>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>>
>>>>>
>>>>> _______________________________________________
>>>>> RsigGeo mailing list
>>>>> [hidden email]
>>>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>>>
>>>>
>>>>
>>>
>>> 
>>> Roger Bivand
>>> Department of Economics, Norwegian School of Economics,
>>> Helleveien 30, N5045 Bergen, Norway.
>>> voice: +47 55 95 93 55; fax +47 55 95 91 00
>>> email: [hidden email]
>>>
>>> _______________________________________________
>>> RsigGeo mailing list
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>
>>
>

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_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo


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: Selfintersection 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: Selfintersection 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 gUnaryUnioning 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: Selfintersection 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.2CAPI1.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, JeanLuc 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 GEOSbased, and returns the correct answer for
>>>> byid=c(TRUE, TRUE) so like rgeos in the same setting. I'm
>>>> reinstalling 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.312 and sp_1.20.
>>>>>>>>
>>>>>>>> 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,
>>>>>>>>
>>>>>>>> JeanLuc Dupouey
>>>>>>>>
>>>>>>>> INRA
>>>>>>>> Forest Ecology & Ecophysiology Unit
>>>>>>>> F54280 Champenoux
>>>>>>>> France
>>>>>>>> mail: [hidden email]
>>>>>>>>
>>>>>>>> Le 23/09/2015 19:39, Roger Bivand a écrit :
>>>>>>>>>
>>>>>>>>> On Wed, 23 Sep 2015, JeanLuc 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,
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>
>>>>>>>> _______________________________________________
>>>>>>>> RsigGeo mailing list
>>>>>>>> [hidden email]
>>>>>>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> _______________________________________________
>>>>>> RsigGeo mailing list
>>>>>> [hidden email]
>>>>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>>>>
>>>>>
>>>>>
>>>>
>>>> 
>>>> Roger Bivand
>>>> Department of Economics, Norwegian School of Economics,
>>>> Helleveien 30, N5045 Bergen, Norway.
>>>> voice: +47 55 95 93 55; fax +47 55 95 91 00
>>>> email: [hidden email]
>>>>
>>>> _______________________________________________
>>>> RsigGeo mailing list
>>>> [hidden email]
>>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>>
>>>
>>
>
>
>
> _______________________________________________
> RsigGeo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/rsiggeo>

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_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo


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: Selfintersection 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.2CAPI1.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, JeanLuc 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 GEOSbased, and returns the correct answer for
>>>> byid=c(TRUE, TRUE) so like rgeos in the same setting. I'm
>>>> reinstalling 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.312 and sp_1.20.
>>>>>>>>
>>>>>>>> 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,
>>>>>>>>
>>>>>>>> JeanLuc Dupouey
>>>>>>>>
>>>>>>>> INRA
>>>>>>>> Forest Ecology & Ecophysiology Unit
>>>>>>>> F54280 Champenoux
>>>>>>>> France
>>>>>>>> mail: [hidden email]
>>>>>>>>
>>>>>>>> Le 23/09/2015 19:39, Roger Bivand a écrit :
>>>>>>>>>
>>>>>>>>> On Wed, 23 Sep 2015, JeanLuc 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,
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>
>>>>>>>> _______________________________________________
>>>>>>>> RsigGeo mailing list
>>>>>>>> [hidden email]
>>>>>>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> _______________________________________________
>>>>>> RsigGeo mailing list
>>>>>> [hidden email]
>>>>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>>>>
>>>>>
>>>>>
>>>>
>>>> 
>>>> Roger Bivand
>>>> Department of Economics, Norwegian School of Economics,
>>>> Helleveien 30, N5045 Bergen, Norway.
>>>> voice: +47 55 95 93 55; fax +47 55 95 91 00
>>>> email: [hidden email]
>>>>
>>>> _______________________________________________
>>>> RsigGeo mailing list
>>>> [hidden email]
>>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>>
>>>
>>
>
> 
> 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>
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo


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 DE9IM 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: Selfintersection 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.2CAPI1.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, JeanLuc 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 GEOSbased, and returns the correct answer for
>>>>> byid=c(TRUE, TRUE) so like rgeos in the same setting. I'm
>>>>> reinstalling 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.312 and sp_1.20.
>>>>>>>>>
>>>>>>>>> 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,
>>>>>>>>>
>>>>>>>>> JeanLuc Dupouey
>>>>>>>>>
>>>>>>>>> INRA
>>>>>>>>> Forest Ecology & Ecophysiology Unit
>>>>>>>>> F54280 Champenoux
>>>>>>>>> France
>>>>>>>>> mail: [hidden email]
>>>>>>>>>
>>>>>>>>> Le 23/09/2015 19:39, Roger Bivand a écrit :
>>>>>>>>>>
>>>>>>>>>> On Wed, 23 Sep 2015, JeanLuc 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,
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>> _______________________________________________
>>>>>>>>> RsigGeo mailing list
>>>>>>>>> [hidden email]
>>>>>>>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> RsigGeo mailing list
>>>>>>> [hidden email]
>>>>>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>> 
>>>>> Roger Bivand
>>>>> Department of Economics, Norwegian School of Economics,
>>>>> Helleveien 30, N5045 Bergen, Norway.
>>>>> voice: +47 55 95 93 55; fax +47 55 95 91 00
>>>>> email: [hidden email]
>>>>>
>>>>> _______________________________________________
>>>>> RsigGeo mailing list
>>>>> [hidden email]
>>>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>>>
>>>>
>>>
>>
>> 
>> 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_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo

Administrator

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: Selfintersection 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: Selfintersection 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 gUnaryUnioning MultiPolygons before intersecting when byid =
> FALSE is not such a bad idea after all?
>
Commited to Rforge 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: Selfintersection 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.2CAPI1.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, JeanLuc 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 GEOSbased, and returns the correct answer for
>>>>> byid=c(TRUE, TRUE) so like rgeos in the same setting. I'm
>>>>> reinstalling 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.312 and sp_1.20.
>>>>>>>>>
>>>>>>>>> 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,
>>>>>>>>>
>>>>>>>>> JeanLuc Dupouey
>>>>>>>>>
>>>>>>>>> INRA
>>>>>>>>> Forest Ecology & Ecophysiology Unit
>>>>>>>>> F54280 Champenoux
>>>>>>>>> France
>>>>>>>>> mail: [hidden email]
>>>>>>>>>
>>>>>>>>> Le 23/09/2015 19:39, Roger Bivand a écrit :
>>>>>>>>>>
>>>>>>>>>> On Wed, 23 Sep 2015, JeanLuc 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,
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>> _______________________________________________
>>>>>>>>> RsigGeo mailing list
>>>>>>>>> [hidden email]
>>>>>>>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> RsigGeo mailing list
>>>>>>> [hidden email]
>>>>>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>> 
>>>>> Roger Bivand
>>>>> Department of Economics, Norwegian School of Economics,
>>>>> Helleveien 30, N5045 Bergen, Norway.
>>>>> voice: +47 55 95 93 55; fax +47 55 95 91 00
>>>>> email: [hidden email]
>>>>>
>>>>> _______________________________________________
>>>>> RsigGeo mailing list
>>>>> [hidden email]
>>>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>>>>
>>>>
>>>
>>
>>
>>
>> _______________________________________________
>> RsigGeo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>
>
>

Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
email: [hidden email]
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N5045 Bergen, Norway

