Overlay between polygons and their intersecting lines

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)

  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.

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)

--
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

--
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