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

[Graham B. McNamara]

Thanks for the referralI. I will start reading up on JTS.

I still don't understand why the example I provided does not run into
problems with this, though. The line and polygon border in my actual
data look similar to the example.





On Thu, Sep 11, 2014 at 11:33 PM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> 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
>