[geos-devel] Intersection Polygon/Polygon curiosity

[Roger Bivand]

Hi,

An oddity has shown up when a toy example is run in R/rgeos and Shapely 
when two polygon geometrycollections are intersected; thread at:

https://stat.ethz.ch/pipermail/r-sig-geo/2015-September/023468.html

With the data as input, GEOSIntersection_r returns the first object 
instead of the intersection. When the coordinates are multiplied by 5, it 
fails with TopologyException: no outgoing dirEdge found ...

Here is the Shapely case (GEOS 3.5.0):

from shapely.wkt import dumps, loads

MyLay = loads("GEOMETRYCOLLECTION (POLYGON ((0.0000000000000000
0.0000000000000000, 0.0000000000000000 1.0000000000000000,
1.0000000000000000 1.0000000000000000, 1.0000000000000000
0.0000000000000000, 0.0000000000000000 0.0000000000000000)), POLYGON
((0.0000000000000000 0.0000000000000000, 1.0000000000000000
0.0000000000000000, 1.0000000000000000 -1.0000000000000000,
0.0000000000000000 -1.0000000000000000, 0.0000000000000000
0.0000000000000000)))")

MyLayl = loads("GEOMETRYCOLLECTION (POLYGON ((0.1000000000000000
0.1000000000000000, 0.1000000000000000 1.1000000000000001,
1.1000000000000001 1.1000000000000001, 1.1000000000000001
0.1000000000000000, 0.1000000000000000 0.1000000000000000)), POLYGON
((0.1000000000000000 0.1000000000000000, 1.1000000000000001
0.1000000000000000, 1.1000000000000001 -0.9000000000000000,
0.1000000000000000 -0.9000000000000000, 0.1000000000000000
0.1000000000000000)))")

merged = MyLay.intersection(MyLayl)

dumps(merged)
'POLYGON ((0.0000000000000000 0.0000000000000000, 0.0000000000000000
1.0000000000000000, 1.0000000000000000 1.0000000000000000,
1.0000000000000000 0.0000000000000000, 1.0000000000000000
-1.0000000000000000, 0.0000000000000000 -1.0000000000000000,
0.0000000000000000 0.0000000000000000))'

merged.equals(MyLay)
True


in R (GEOS 3.5.0):

library(sp)
library(rgeos)
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))

Pols1=Polygons(list(Polygon(Pol1)),"Pols1")
Pols2=Polygons(list(Polygon(Pol2)),"Pols2")
MyLay=SpatialPolygons(list(Pols1,Pols2))


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

inter=gIntersection(MyLay,MyLayl)
gEquals(inter, MyLay)

plot(MyLay)
plot(MyLayl,add=TRUE)
plot(inter,col="red",add=TRUE)
writeWKT(MyLay)
writeWKT(MyLayl)

There is a report on the thread of PostGIS not showing the issue, but we 
are checking that the same input data are being used.

Grateful for any input correcting my misunderstanding ...

Roger

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