[R-sig-Geo] gIntersection returns error "TopologyException: no outgoing dirEdge found at"
Edzer Pebesma
edzer.pebesma at uni-muenster.de
Fri Sep 25 22:26:44 CEST 2015
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?
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
>
--
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/49b1dcc0/attachment.bin>
More information about the R-sig-Geo
mailing list