[R-sig-Geo] gIntersection with drop_lower_td=T drops entire intersection

[Roger Bivand]

Pius,

There was a bug in the C code which caused the deletion of the Polygons 
object rather than the offending point in the geometry collection. The 
current source of rgeos on R-forge should fix the problem - I submitted 
source tarball the source tarball to R Winbuilder, and a Windows binary 
build of the development package is available on (if you use Windows, R 
3.2.2):

http://win-builder.r-project.org/WWjSD9bDw8xx

Thanks for your report of this problem, your example code pinpointed the 
issue, and is now among the examples on the help page.

Best wishes,

Roger

On Fri, 23 Oct 2015, Roger Bivand wrote:

> I didn't see this when you posted, sorry. I can replicate the problem - 
> thanks for a useful report. I'll try to see where drop_lower_td= is being too 
> aggresive.
>
> I'd be very obliged if others (not necessarily the poster!) would try to 
> answer questions like this, the code is open, my time is limited, and 
> sometimes my day-job has priority ...
>
> Roger
>
> On Thu, 1 Oct 2015, Pius Korner wrote:
>
>> Dear all
>> 
>> 
>> 
>> I intersect a SpatialPolygons-Object with two polygons 1a and 1b with a
>> polygon 2. Polygon 1a intersects with polygon 2 in one area but also only
>> touches polygon 2 at another position.
>> 
>> Here a schematic case for illustration:
>> 
>> 
>> 
>> poly1 <-
>> SpatialPolygons(list(Polygons(list(Polygon(coords=matrix(c(0,0,2,2,0,1,1,0),
>> ncol=2,byrow=F))),ID=c("a")),
>> 
>> 
>> Polygons(list(Polygon(coords=matrix(c(0,0,2,2,2,3,3,2),ncol=2,byrow=F))),ID=
>> c("b"))))
>> 
>> poly2 <-
>> SpatialPolygons(list(Polygons(list(Polygon(coords=matrix(c(0,0,2,2,1,1,1,3,3
>> ,0,0,2),ncol=2,byrow=F))),ID=c("c"))))
>> 
>> plot(poly1,border="orange")
>> 
>> plot(poly2,border="blue",add=T,lty=2,density=8,angle=30,col="blue")
>> 
>> 
>> 
>> I'm interested in the area of overlap of polygons 1a and 1b with polygon 2.
>> 
>> gIntersection(poly1,poly2,byid=T)
>> 
>> does not work because of the different subgeometries that result from the
>> intersection of 1a with 2; they cannot be collected.
>> 
>> gIntersection(poly1,poly2,byid=T,drop_lower_td=T)
>> 
>> drops the entire intersection between 1a with 2, hence I get no value for
>> 1a, even though it does overlap with polygon 2.
>> 
>> 
>> 
>> Is there a way so that only the Point-subgeometries (with area=0) are
>> dropped, but the Polygon-subgeometry is retained and returned for the
>> intersection of 1a and 2?
>> 
>> 
>> 
>> my temporary workaround was to very slightly buffer poly1 (width=-0.0001),
>> to prevent the point-subgeometry.
>> 
>> I could create a loop and intersect each polygon in poly1 separately, but
>> using byid=T would be nicer (faster..)
>> 
>> 
>> 
>> many thanks
>> 
>> Pius
>> 
>>
>> 	[[alternative HTML version deleted]]
>> 
>> _______________________________________________
>> 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