On 28/05/16 00:21, Vijay Lulla wrote:
Edzer,
I'm having trouble replicating what you state. Below is from my R session.
I don't see what it is you have trouble with (except lacking documentation)
### Start R session interaction
over(Inter2[1], PolyD2, byid=T)
# The green (1st) segment intersects with the black (1st) polygon => good
over(Inter2[2], PolyD2, byid=T) # Still not ok!
# 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!
over(Inter2[2], PolyD2, byid=T, returnList = TRUE, minDimension = 1) # OK!
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.
Sorry, my mistake. You'll find details in vignette("over")
over() finds intersections without ordering them, so if a line overlaps
polygon A but touches polygon B, it might as well return polygon B as
your match, where you were hoping for polygon A.
With minDimension = 1 the dimensionality of the intersection (read the
help of rgeos::gRelate) needs to be at least one, so touching (a point:
0-dimensional) is ignored, and polygon A would be returned.
minDimension is now 0 by default; other values might make more sense,
depending on the geometries at hand. A problem is that this increases
the computational burden quite a bit (requiring gRelate instead of
gIntersects).
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)