[R-sig-Geo] gIntersection returns error "TopologyException: no outgoing dirEdge found at"
Edzer Pebesma
edzer.pebesma at uni-muenster.de
Fri Sep 25 21:55:20 CEST 2015
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
>>>
>>
>
--
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
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 490 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20150925/2d036967/attachment.bin>
More information about the R-sig-Geo
mailing list