[R-sig-Geo] Overlay between polygons and their intersecting lines

[Vijay Lulla]

Edzer,
I'm having trouble replicating what you state.  Below is from my R session.

### Start R session interaction

> over(Inter2[1], PolyD2, byid=T)
  Attribute
1    1black
> # The green (1st) segment intersects with the black (1st) polygon => good
> over(Inter2[2], PolyD2, byid=T)  # Still not ok!
  Attribute
1    1black
> # The blue(2nd) segment intersects with the red(2nd) polygon => NOW GOOD !
# The blue(2nd) segment intersects with the red(2nd) polygon => NOW GOOD !
> over(Inter2[1], PolyD2, byid=T, returnList = TRUE, minDimension = 1) # OK!
[[1]]
  Attribute
1    1black

> over(Inter2[2], PolyD2, byid=T, returnList = TRUE, minDimension = 1) # OK!
[[1]]
  Attribute
2      2red

> over(Inter2[1], PolyD2, byid=T, returnList=TRUE)
>
>
> over(Inter2[1], PolyD2, byid=T, returnList=TRUE)
[[1]]
  Attribute
1    1black
2      2red

> over(Inter2[2], PolyD2, byid=T, returnList=TRUE)
[[1]]
  Attribute
1    1black
2      2red

### End R session

By the way, what is the minDimension argument in `over` function?  Why
do we need it?  Also, where can I read about it?  I cannot find it in
?over.

Thanks,
Vijay.

On Fri, May 27, 2016 at 5:33 PM, Edzer Pebesma
<edzer.pebesma at uni-muenster.de> wrote:
> Thanks for the clear example. I can reproduce this, but don't know how
> to solve it. If I create a similar example on the unit square, things
> work fine. The different outputs of geos::gRelate() indicate the problem
> (script below):
>
>> rgeos::gRelate(Inter, PolyD, byid = TRUE) # your example:
>   1 1         1 2
> 1 "1FF00F212" "1010F0212"
> 2 "FF1FF0212" "101F00212"
>
>> rgeos::gRelate(Inter2, PolyD2, byid = TRUE) # simplified numbers:
>   1 1         1 2
> 1 "1FFF0F212" "FF1F00212"
> 2 "FF1F00212" "1FFF0F212"
>
> Is it possible that this relates to the scale factor of GEOS?
>
> # script:
> library(sp)
> matrix1 <- cbind(x= c(0, 1, 1, 0), y= c(2, 2, 0, 0)) # Simplified!
> matrix2 <- cbind(x= c(1, 2, 2, 1), y= c(2, 2, 0, 0))
>
> par(mfrow = c(2,1))
>
> P1 <- Polygon (matrix1)
> P2 <- Polygon (matrix2)
> Ps1 <- Polygons (list(P1), 1 )
> Ps2 <- Polygons (list(P2), 2 )
> Poly <- SpatialPolygons (list (Ps1, Ps2))
> data <- data.frame( Attribute = c("1black", "2red") )
> PolyD2 <- SpatialPolygonsDataFrame(Poly , data = data )
> matrix3 <- cbind(x= c(0, 2), y= c(0.5, 0.5))
> L <- Line(matrix3)
> Ls <- Lines(list(L), "1")
> Line <- SpatialLines( list(Ls))
> Inter2 <- rgeos::gIntersection(Line,PolyD2, byid=T)
> plot(PolyD2, col = PolyD2$Attribute, axes =TRUE)
> lines( Inter2[1], col = "green", lwd = 3 )
> lines( Inter2[2], col = "blue", lwd = 3 )
> over(Inter2[1], PolyD2, byid=T)
> # The green (1st) segment intersects with the black (1st) polygon => good
> over(Inter2[2], PolyD2, byid=T)
> # The blue(2nd) segment intersects with the red(2nd) polygon => NOW GOOD !
> over(Inter2[1], PolyD2, byid=T, returnList = TRUE, minDimension = 1) # OK!
> over(Inter2[2], PolyD2, byid=T, returnList = TRUE, minDimension = 1) # OK!
>
> matrix1 <- cbind(x= c(250300, 250451, 250494, 250300), y= c(104030,
> 104030, 103733, 103733))
> matrix2 <- cbind(x= c(250451, 250666, 250666, 250494), y= c(104030,
> 104030, 103733, 103733))
> P1 <- Polygon (matrix1)
> P2 <- Polygon (matrix2)
> Ps1 <- Polygons (list(P1), 1 )
> Ps2 <- Polygons (list(P2), 2 )
> Poly <- SpatialPolygons (list (Ps1, Ps2))
> data <- data.frame( Attribute = c("1black", "2red") )
> PolyD <- SpatialPolygonsDataFrame(Poly , data = data )
> plot(PolyD, col = PolyD$Attribute, axes =TRUE)
> matrix3 <- cbind(x= c(250300, 250666),
> y= c(103900, 103900))
> L <- Line(matrix3)
> Ls <- Lines(list(L), "1")
> Line <- SpatialLines( list(Ls))
> Inter <- rgeos::gIntersection(Line,PolyD, byid=T)
> lines( Inter[1], col = "green", lwd = 3 )
> lines( Inter[2], col = "blue", lwd = 3 )
> over(Inter[1], PolyD, byid=T)
> # The green (1st) segment intersects with the black (1st) polygon => good
> over(Inter[2], PolyD, byid=T)
> # The blue(2nd) segment intersects with the red(2nd) polygon => NOT GOOD !
> over(Inter[1], PolyD, byid=T, returnList = TRUE) # OK
> over(Inter[2], PolyD, byid=T, returnList = TRUE)
>
> # illustrate the problem:
> rgeos::gRelate(Inter, PolyD, byid = TRUE)
> rgeos::gRelate(Inter2, PolyD2, byid = TRUE)
>
>
> There's also an in-line comment below:
>
>
> On 25/05/16 15:53, Gwennaël Bataille wrote:
>> Dear all,
>>
>> I can't find a solution for the following problem:
>>
>> When I first intersect a line with 2 polygons (splitting it into 2
>> segments) and then use an overlay to get for each segment the attribute
>> of the overlapping polygon, I sometimes get too answers (i.e. a small
>> point overlapping one polygon, the rest of the segment overlapping another).
>>
>> The functions I use for this are:
>>
>> rgeos::gIntersection
>>
>> and sp::over
>>
>> Do you have any idea how I could get the attribute of the polygon the
>> segments are "mostly overlapping"?
>>
>> Below is a reproductible example.
>>
>> Thank you very much in advance for any hint on this.
>>
>> Best regards,
>>
>> Gwenna�l
>>
>>
>>
>> Reproductible example:
>>
>> matrix1 <- cbind(x= c(250300, 250451, 250494, 250300),
>>
>> y= c(104030, 104030, 103733, 103733))
>>
>> matrix2 <- cbind(x= c(250451, 250666, 250666, 250494),
>>
>> y= c(104030, 104030, 103733, 103733))
>>
>> P1 <- Polygon (matrix1)
>>
>> P2 <- Polygon (matrix2)
>>
>> Ps1 <- Polygons (list(P1), 1 )
>>
>> Ps2 <- Polygons (list(P2), 2 )
>>
>> Poly <- SpatialPolygons (list (Ps1, Ps2), proj4string=CRS("+proj=longlat
>> +datum=WGS84"))
>
> this proj4string doesn't make sense for such coordinates.
>
>>
>> data <- data.frame( Attribute = c("1black", "2red") )
>>
>> PolyD <- SpatialPolygonsDataFrame(Poly , data = data )
>>
>> plot(PolyD, col = PolyD$Attribute)
>>
>> matrix3 <- cbind(x= c(250300, 250666),
>>
>> y= c(103900, 103900))
>>
>> L <- Line(matrix3)
>>
>> Ls <- Lines(list(L), "1")
>>
>> Line <- SpatialLines( list(Ls), proj4string = CRS("+proj=longlat
>> +datum=WGS84") )
>>
>> Inter <- rgeos::gIntersection(Line,PolyD, byid=T)
>>
>> lines( Inter[1], col = "green", lwd = 3 )
>>
>> lines( Inter[2], col = "blue", lwd = 3 )
>>
>> over(Inter[1], PolyD, byid=T)
>>
>> 'Attribute
>>
>> 11black'
>>
>> # The green (1st) segment intersects with the black (1st) polygon => good
>>
>> over(Inter[2], PolyD, byid=T)
>>
>> 'Attribute
>>
>> 11black'
>>
>> # The blue(2nd) segment intersects with the red(2nd) polygon => NOT GOOD !
>>
>> over(Inter[1], PolyD, byid=T, returnList = TRUE) # OK
>>
>> over(Inter[2], PolyD, byid=T, returnList = TRUE)
>>
>> '[[1]]
>>
>> Attribute
>>
>> 11black
>>
>> 22red'
>>
>> # The blue(2nd) segment appears to intersect the two polygons (I guess
>> one point in common with the first polygon; the rest of the segment
>> overlapping the other polygon)
>>
>>
>>
>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> --
> Edzer Pebesma
> Institute for Geoinformatics  (ifgi),  University of Münster
> Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
> Journal of Statistical Software:   http://www.jstatsoft.org/
> Computers & Geosciences:   http://elsevier.com/locate/cageo/
> Spatial Statistics Society http://www.spatialstatistics.info
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo