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