[R-sig-Geo] gIntersection returns error "TopologyException: no outgoing dirEdge found at"
Vijay Lulla
vijaylulla at gmail.com
Fri Sep 25 19:28:02 CEST 2015
Seems to be working in PostGIS! Below is from my PostGIS session:
gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 10, 10 10, 10 0, 0
0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
st_intersects
---------------
t
(1 row)
gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 10, 10 10, 10 0, 0
0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
st_intersection
------------------------------------------------------------------------------------
0102000000020000000000000000002440000000000000000000000000000000000000000000000000
(1 row)
gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 1, 1 1, 1 0, 0
0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
st_intersection
------------------------------------------------------------------------------------
010200000002000000000000000000F03F000000000000000000000000000000000000000000000000
(1 row)
gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 1, 1 1, 1 0, 0
0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
st_intersects
---------------
t
(1 row)
gisdb=# select postgis_version();
postgis_version
---------------------------------------
2.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1
(1 row)
gisdb=# select postgis_full_version();
postgis_full_vers
ion
------------------------------------------------------------------------------------------
----------------------------------------------------------------------------
POSTGIS="2.1.7 r13414" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel.
4.8.0, 6 March 2012" GDAL
="GDAL 1.11.1, released 2014/09/24" LIBXML="2.7.8" LIBJSON="UNKNOWN" RASTER
(1 row)
On Fri, Sep 25, 2015 at 1:02 PM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> On Fri, 25 Sep 2015, Edzer Pebesma wrote:
>
>>
>>
>> On 25/09/15 14:36, Roger Bivand wrote:
>>>
>>> On Fri, 25 Sep 2015, Roger Bivand wrote:
>>>
>>>> On Fri, 25 Sep 2015, Jean-Luc Dupouey wrote:
>>>>
>>> ...
>>>>>
>>>>>
>>>>> library(sp)
>>>>> library(rgeos)
>>>>>
>>>>> #build a first polygon, MyLay, composed of two adjacent squares, size
>>>>> 1x1:
>>>>>
>>>>> Pol1=rbind(c(0,0),c(0,1),c(1,1),c(1,0),c(0,0))
>>>>> Pol2=rbind(c(0,0),c(1,0),c(1,-1),c(0,-1),c(0,0))
>>>
>>>
>>> Changing Pol1 and Pol2 by * 5, and the increment below to 0.5 gives the
>>> earlier:
>>>
>>>> inter=gIntersection(MyLay,MyLayl)
>>>
>>> Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td,
>>> "rgeos_intersection") :
>>> TopologyException: no outgoing dirEdge found at 0.5 0
>>>
>>> so this is another rendering of the original edge case in GEOS in this
>>> thread. The practical resolution is to run gUnaryUnion() on the objects
>>> when byid=FALSE and the required output is a single (possibly multipart)
>>> geometry containing the whole intersection.
>>>
>>> Should we consider using gUnaryUnion() by default if byid for the object
>>> is FALSE, but permit users to choose not to do this?
>>
>>
>> Has someone tried to reproduce this with GEOS, without R, e.g. by a
>> PostGIS query?
>
>
> Do you have a running PostGIS instance - I don't. I did try v.overlay in
> GRASS, but it isn't GEOS-based, and returns the correct answer for
> byid=c(TRUE, TRUE) so like rgeos in the same setting. I'm re-installing GEOS
> with Python, but someone else should check (my GEOS 3.5.0).
>
> Roger
>
>
>>
>>>
>>> Roger
>>>
>>>>>
>>>>> Pols1=Polygons(list(Polygon(Pol1)),"Pols1")
>>>>> Pols2=Polygons(list(Polygon(Pol2)),"Pols2")
>>>>> MyLay=SpatialPolygons(list(Pols1,Pols2))
>>>>>
>>>>> #build a second polygon, MyLayl, composed of the same two squares
>>>>> lagged by 0.1 in x and y directions:
>>>>>
>>>>> Pol1l=Pol1+0.1
>>>>> Pol2l=Pol2+0.1
>>>>>
>>>>> Pols1l=Polygons(list(Polygon(Pol1l)),"Pols1l")
>>>>> Pols2l=Polygons(list(Polygon(Pol2l)),"Pols2l")
>>>>> MyLayl=SpatialPolygons(list(Pols1l,Pols2l))
>>>>>
>>>>> #view the resulting spatial layers:
>>>>>
>>>>> plot(MyLay)
>>>>> plot(MyLayl,add=TRUE)
>>>>>
>>>>> #"successful" intersection:
>>>>>
>>>>> inter=gIntersection(MyLay,MyLayl)
>>>>
>>>>
>>>> You did not follow my advice to put the objects into gUnaryUnion()
>>>> first, to remove the problem. It isn't obvious where this problem is
>>>> coming from, but most likely requires debugging in GEOS itself.
>>>>
>>>>>
>>>>> #but false result: the intersection gives the same contour as MyLay!
>>>>>
>>>>> plot(MyLay)
>>>>> plot(MyLayl,add=TRUE)
>>>>> plot(inter,col="red",add=TRUE)
>>>>>
>>>>> I use R version 3.2.2, rgeos_0.3-12 and sp_1.2-0.
>>>>>
>>>>> If such an error occurs when crossing larger spatial layers, it will
>>>>> remain undetectable.
>>>>>
>>>>> If I am not wrong, it implies that this software is not very
>>>>> reliable. Especially when considering that byid=FALSE is the default
>>>>> option.
>>>>>
>>>>> gIntersection gives a correct result using the option byid=TRUE.
>>>>
>>>>
>>>> This was not what you said you wanted, which was the (single) object
>>>> for which MyLay intersected MyLayl.
>>>>
>>>>> The intersect function of the raster package also gives a correct
>>>>> intersection. Do you think this latter function could be a better
>>>>> option, in general?
>>>>>
>>>>
>>>> No, not at all. If you inspect it, you'll see that it simply tries to
>>>> "help" users by trying to guess which "data" slot elements might be
>>>> copied across. It also calls gIntersect() in dense mode, which will
>>>> choke on intersections between objects with many geometries. It runs
>>>> gIntersection(..., byid=TRUE), so adds little to just doing that.
>>>>
>>>> See:
>>>>
>>>> library(rgdal)
>>>> getMethod("intersect", c("SpatialPolygons", "SpatialPolygons"))
>>>>
>>>> Hope this clarifies,
>>>>
>>>> Roger
>>>>
>>>>> Thanks in advance,
>>>>>
>>>>> Jean-Luc Dupouey
>>>>>
>>>>> INRA
>>>>> Forest Ecology & Ecophysiology Unit
>>>>> F-54280 Champenoux
>>>>> France
>>>>> mail: dupouey at nancy.inra.fr
>>>>>
>>>>> Le 23/09/2015 19:39, Roger Bivand a écrit :
>>>>>>
>>>>>> On Wed, 23 Sep 2015, Jean-Luc Dupouey wrote:
>>>>>>
>>>>>>> Why gIntersection returns the ""TopologyException: no outgoing
>>>>>>> dirEdge found at" error in the following very simple case:
>>>>>>>
>>>>>>>> #construct a spatial layer with two adjacent squares, each of 10 x
>>>>>>>> 10
>>>>>>>
>>>>>>> units:
>>>>>>>>
>>>>>>>>
>>>>>>>> Pol1=rbind(c(0,0),c(0,10),c(10,10),c(10,0),c(0,0))
>>>>>>>> Pol2=rbind(c(0,0),c(10,0),c(10,-10),c(0,-10),c(0,0))
>>>>>>>>
>>>>>>>> Pols1=Polygons(list(Polygon(Pol1)),"Pols1")
>>>>>>>> Pols2=Polygons(list(Polygon(Pol2)),"Pols2")
>>>>>>>> MyLay=SpatialPolygons(list(Pols1,Pols2))
>>>>>>>>
>>>>>>>> #construct the same layer, but lagged by 0.5 in x and y directions
>>>>>>>>
>>>>>>>> Pol1l=Pol1+0.5
>>>>>>>> Pol2l=Pol2+0.5
>>>>>>>>
>>>>>>>> Pols1l=Polygons(list(Polygon(Pol1l)),"Pols1l")
>>>>>>>> Pols2l=Polygons(list(Polygon(Pol2l)),"Pols2l")
>>>>>>>> MyLayl=SpatialPolygons(list(Pols1l,Pols2l))
>>>>>>>>
>>>>>>>> #view the resulting spatial layers:
>>>>>>>>
>>>>>>>> plot(MyLay)
>>>>>>>> plot(MyLayl,add=TRUE)
>>>>>>>>
>>>>>>>> #try to intersect:
>>>>>>>>
>>>>>>>> inter=gIntersection(MyLay,MyLayl)
>>>>>>>
>>>>>>> Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td,
>>>>>>> "rgeos_intersection") :
>>>>>>> TopologyException: no outgoing dirEdge found at 0.5 0
>>>>>>
>>>>>>
>>>>>> The error message is coming from GEOS - you are welcome to investigate
>>>>>> further. If you use gUnaryUnion() on the objects first, there is no
>>>>>> error:
>>>>>>
>>>>>> library(rgeos)
>>>>>> Pol1=rbind(c(0,0),c(0,10),c(10,10),c(10,0),c(0,0))
>>>>>> Pol2=rbind(c(0,0),c(10,0),c(10,-10),c(0,-10),c(0,0))
>>>>>> library(sp)
>>>>>> Pols1=Polygons(list(Polygon(Pol1)),"Pols1")
>>>>>> Pols2=Polygons(list(Polygon(Pol2)),"Pols2")
>>>>>> MyLay=SpatialPolygons(list(Pols1,Pols2))
>>>>>> Pol1l=Pol1+0.5
>>>>>> Pol2l=Pol2+0.5
>>>>>> Pols1l=Polygons(list(Polygon(Pol1l)),"Pols1l")
>>>>>> Pols2l=Polygons(list(Polygon(Pol2l)),"Pols2l")
>>>>>> MyLayl=SpatialPolygons(list(Pols1l,Pols2l))
>>>>>> plot(MyLay)
>>>>>> plot(MyLayl,add=TRUE)
>>>>>> points(x=0.5, y=0)
>>>>>> inter=gIntersection(gUnaryUnion(MyLay), gUnaryUnion(MyLayl))
>>>>>> plot(inter, add=TRUE, border="green")
>>>>>>
>>>>>> It is a GEOS issue, not rgeos.
>>>>>>
>>>>>>>
>>>>>>> It works with the option byid=TRUE.
>>>>>>>
>>>>>>> But my question is why it does not work without? Is this behaviour
>>>>>>> predictable?
>>>>>>
>>>>>>
>>>>>> This is unknown. I'll try again on GEOS 3.5.0 and see if it is still
>>>>>> there: yes, unfortunately it is still there.
>>>>>>
>>>>>> Roger
>>>>>>
>>>>>>>
>>>>>>> I went through some previous related posts to the list, but could not
>>>>>>> find the answer.
>>>>>>>
>>>>>>> Thanks,
>>>>>>>
>>>>>>>
>>>>>>
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-Geo mailing list
>>>>> R-sig-Geo at r-project.org
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>
>>>>
>>>>
>>>
>>>
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>>
>
> --
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; fax +47 55 95 91 00
> e-mail: Roger.Bivand at nhh.no
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
More information about the R-sig-Geo
mailing list