Skip to content

gIntersection with drop_lower_td=T drops entire intersection

3 messages · Pius Korner, Roger Bivand

Original
Pius Korner
# 
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
21 days later
Roger Bivand
# 
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:

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    
1 day later

    

      
      

        


  

        

      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: