[R-sig-Geo] gIntersection gives error (may not contain geometry collections)
Roger.Bivand at nhh.no
Thu Sep 11 23:33:33 CEST 2014
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,
> 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?
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
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