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

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,

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