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

[Edzer Pebesma]

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

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 490 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20160527/1a5abaf8/attachment.bin>