[R-sig-Geo] gIntersection returns error "TopologyException: no outgoing dirEdge found at"
Roger Bivand
Roger.Bivand at nhh.no
Sat Sep 26 19:29:09 CEST 2015
On Fri, 25 Sep 2015, Edzer Pebesma wrote:
> Trying to reduce the problem, it seems PostGIS thinks that polygons are
> invalid for which R believes they are valid:
>
> postgis=# select ST_IsValid('MULTIPOLYGON(
> postgis'# ((0 0, 0 10, 10 10, 10 0, 0 0)),
> postgis'# ((0 0, 10 0, 10 -10, 0 -10, 0 0)))'::geometry);
> NOTICE: Self-intersection at or near point 0 0
> st_isvalid
> ------------
> f
> (1 row)
>
> postgis=# select ST_IsValid('POLYGON(
> postgis'# (0 0, 0 10, 10 10, 10 0, 0 0),
> postgis'# (0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
> NOTICE: Self-intersection at or near point 0 0
> st_isvalid
> ------------
> f
> (1 row)
>
> but
>
>> library(rgeos)
>> Pol1=rbind(c(0,0),c(0,10),c(10,10),c(10,0),c(0,0))
>> Pol2=rbind(c(0,0),c(10,0),c(10,-10),c(0,-10),c(0,0))
>>
>> Pols1=Polygons(list(Polygon(Pol1)),"Pols1")
>> Pols2=Polygons(list(Polygon(Pol2)),"Pols2")
>> MyLay=SpatialPolygons(list(Pols1,Pols2))
>> gIsValid(MyLay)
> [1] TRUE
>
> Maybe gUnaryUnion-ing MultiPolygons before intersecting when byid =
> FALSE is not such a bad idea after all?
>
Commited to R-forge as revision 505, seems to address the case, please
check. The fix doesn't branch on geometry types, could those affected
please check whether this leads to regressions?
Roger
>
> On 25/09/15 21:55, Edzer Pebesma wrote:
>> Thanks Vijay; FYI, following OP's example, I see in PostGIS:
>>
>> SELECT ST_AsText(ST_Intersection(
>> 'MULTIPOLYGON(
>> ((0 0, 0 10, 10 10, 10 0, 0 0)),
>> ((0 0, 10 0, 10 -10, 0 -10, 0 0)))'::geometry,
>> 'MULTIPOLYGON(
>> ((0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5)),
>> ((0.5 0.5, 10.5 0.5, 10.5 -10.5, 0.5 -10.5, 0.5 0.5)))'::geometry
>> ));
>>
>> giving:
>>
>> ERROR: Error performing intersection: TopologyException: Input geom 0
>> is invalid: Self-intersection at or near point 0 0 at 0 0
>>
>> but
>>
>> SELECT ST_AsText(ST_Intersection(
>> 'MULTIPOLYGON(
>> ((0 0, 0 10, 10 10, 10 0, 0 0)),
>> ((0 -1, 10 -1, 10 -10, 0 -10, 0 -1)))'::geometry,
>> 'MULTIPOLYGON(
>> ((0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5)),
>> ((0.5 0.0, 10.5 0.0, 10.5 -10.5, 0.5 -10.5, 0.5 0.0)))'::geometry
>> ));
>>
>> giving no error; when I try
>>
>> SELECT ST_AsText(ST_Intersection(
>> 'POLYGON(
>> (0 0, 0 10, 10 10, 10 0, 0 0),
>> (0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry,
>> 'POLYGON(
>> (0.5 0.5, 0.5 10.5, 10.5 10.5, 10.5 0.5, 0.5 0.5),
>> (0.5 0.5, 10.5 0.5, 10.5 -10.5, 0.5 -10.5, 0.5 0.5))'::geometry
>> ));
>>
>> it gives the same error as the first example.
>>
>>
>>
>> On 25/09/15 19:44, Roger Bivand wrote:
>>> On Fri, 25 Sep 2015, Vijay Lulla wrote:
>>>
>>>> Seems to be working in PostGIS! Below is from my PostGIS session:
>>>>
>>>> gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 10, 10 10, 10 0, 0
>>>> 0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
>>>> st_intersects
>>>> ---------------
>>>> t
>>>> (1 row)
>>>>
>>>>
>>>> gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 10, 10 10, 10 0, 0
>>>> 0))'::geometry, 'POLYGON((0 0, 10 0, 10 -10, 0 -10, 0 0))'::geometry);
>>>> st_intersection
>>>> ------------------------------------------------------------------------------------
>>>>
>>>> 0102000000020000000000000000002440000000000000000000000000000000000000000000000000
>>>>
>>>> (1 row)
>>>>
>>>> gisdb=# SELECT ST_Intersection('POLYGON((0 0, 0 1, 1 1, 1 0, 0
>>>> 0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
>>>> st_intersection
>>>> ------------------------------------------------------------------------------------
>>>>
>>>> 010200000002000000000000000000F03F000000000000000000000000000000000000000000000000
>>>>
>>>> (1 row)
>>>
>>> Thanks for this, but are these the same - that is each geometry
>>> collection consists of two touching squares, but the second is offset
>>> =.1 North Eastwards? In any case, very grateful that you took a look!
>>>
>>> Roger
>>>
>>>>
>>>>
>>>> gisdb=# SELECT ST_Intersects('POLYGON((0 0, 0 1, 1 1, 1 0, 0
>>>> 0))'::geometry, 'POLYGON((0 0, 1 0, 1 -1, 0 -1, 0 0))'::geometry);
>>>> st_intersects
>>>> ---------------
>>>> t
>>>> (1 row)
>>>>
>>>> gisdb=# select postgis_version();
>>>> postgis_version
>>>> ---------------------------------------
>>>> 2.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1
>>>> (1 row)
>>>>
>>>>
>>>> gisdb=# select postgis_full_version();
>>>>
>>>> postgis_full_vers
>>>> ion
>>>> ------------------------------------------------------------------------------------------
>>>>
>>>> ----------------------------------------------------------------------------
>>>>
>>>> POSTGIS="2.1.7 r13414" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel.
>>>> 4.8.0, 6 March 2012" GDAL
>>>> ="GDAL 1.11.1, released 2014/09/24" LIBXML="2.7.8" LIBJSON="UNKNOWN"
>>>> RASTER
>>>> (1 row)
>>>>
>>>> On Fri, Sep 25, 2015 at 1:02 PM, Roger Bivand <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
>>>>>
>>>>
>>>
>>
>>
>>
>> _______________________________________________
>> 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