[R-sig-Geo] gIntersection gives error (may not contain geometry collections)

[Roger Bivand]

On Thu, 11 Sep 2014, Graham Bond McNamara wrote:

> I am trying to use gIntersection like so:
>
> gIntersection(mypolys.spdf, mylines.sldf, byid=c(TRUE,TRUE))
>
> This intersection between the polygons in the SpatialPolygonsDataFrame
> and the lines in the SpatialLinesDataFrame should result in a new
> SpatialLines object. I have used this before and this works.
>
> However, on my new dataset, I get the following error:
>
> output subgeometry 3195, row.name: 24 3759
> subsubgeometry 0: Point
> subsubgeometry 1: LineString
> Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_not_poly,
> "rgeos_intersection") :
>  Geometry collections may not contain other geometry collections
>
> I have inspected mypolys.spdf and my lines.sldf, and they are not
> geometry collections, but just the spdf and sldf. So I don't know
> where the "Point" and "LineString" are coming from (what is a
> LineString anyway?).

Briefly, an input Polygons object and an input Lines object intersect at a 
point and a line, so returning a collection of objects of different types. 
Almost certainly the intention is to drop the point (a line of no length), 
as is already optionally done when a polygon/polygon intersection returns 
polygon and non-polygon (point or line) objects - the drop_not_poly= 
argument. Your case suggests that this should be checked again, to change 
the present argument to drop points for lines and points/lines for 
polygons. The details of the meanings of operations are found in the JTS 
references in the package.

Hope this clarifies,

Roger

>
> The culprit seems to be the combination of  poly 24 and line 3759.
> There is nothing wrong with it per se, but it's a line that exactly
> overlaps the polygon border. I thought perhaps that was the problem?
>
> I am not allowed to share my actual data, but it looks like this:
>
> Sr1 = Polygon(cbind(c(1,2,2,4,4,1,1),c(1,1,3,3,5,5,1)))
> Sr2 = Polygon(cbind(c(2,4,4,2,2),c(1,1,3,3,1)))
> Srs1 = Polygons(list(Sr1), "1")
> Srs2 = Polygons(list(Sr2), "2")
> SpP = SpatialPolygons(list(Srs1,Srs2), 1:2)
> plot(SpP) # example neighborhoods
>
> l1 = cbind(c(2,2),c(2,4))
> Sl1 = Line(l1)
> S1 = Lines(list(Sl1), ID="1")
> Sl = SpatialLines(list(S1))
> plot(Sl, add=TRUE, col="green", lwd=3)
>
> gIntersection works:
>
> gInt <- gIntersection(SpP, Sl, byid=c(TRUE, TRUE))
>
> The line hugging the border intersects both poly 1 and poly 2, so this
> line is split up and reproduced twice.
>
> In summary: gIntersection works even when a line is exactly on the
> poly's border. This means I am no closer to understanding my error
> message. Can anyone help?
>
> Cheers,
> Graham
>
> _______________________________________________
> 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