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

[Gwennaël Bataille]

Thank you very much for your answers!

rgeos::gRelate(Inter[1], PolyD, byid = TRUE)
# From what I undestand of the "Dimensionally Extended nine-Intersection 
Model" ( https://en.wikipedia.org/wiki/DE-9IM ),
# It looks like the green segment is in fact a bit shorter and the blue 
one starts a bit before the separation between the two polygons
# Then it makes sense (e.g. the intersection between the boundary of the 
green segment and the interior of the black polygon is the end point,
# there is no intersection between the boundaries of the green segment 
and the red polygon, ...)

"Is it possible that this relates to the scale factor of GEOS?"
# It seems so; thanks for the hint!
# Changing the scale before gIntersection() changes the outcome.
# Alternatively, we can change it before over(). Strangely, I only found 
one "good scale" for this, though I don't know how to generalize the 
result (which scale to apply when).

# Only "scale = 0.01" leads to the good result:
setScale( scale=0.01 )
#Inter1b <- gDifference( Inter[1], gBoundary(Inter[1]) )
#Inter2b <- gDifference( Inter[2], gBoundary(Inter[2]) )
over(Inter[1], PolyD, byid=T, returnList = TRUE, minDimension = 1) # OK
[[1]]
   Attribute
1    1black
over(Inter[2], PolyD, byid=T, returnList = TRUE, minDimension = 1) # OK
[[1]]
   Attribute
2      2red



# Wrong scales for the over() output:

setScale( scale=0.01 / 10 )
over(Inter[1], PolyD, byid=T, returnList = TRUE, minDimension = 1) # nothing
over(Inter[2], PolyD, byid=T, returnList = TRUE, minDimension = 1) # 2 
overlaps

setScale( scale=1 )
over(Inter[1], PolyD, byid=T, returnList = TRUE, minDimension = 1) # 2 
overlaps
over(Inter[2], PolyD, byid=T, returnList = TRUE, minDimension = 1) # OK

setScale( scale=10 )
over(Inter[1], PolyD, byid=T, returnList = TRUE, minDimension = 1) # OK
over(Inter[2], PolyD, byid=T, returnList = TRUE, minDimension = 1) # 2 
overlaps



All the best,


Gwennael



Le 28/05/2016 12:15, Edzer Pebesma a écrit :
>
> 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)
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> _______________________________________________
>>>> 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
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

-- 
Gwennaël BATAILLE, PhD student - Teaching assistant

Earth and Life Institute
Université Catholique de Louvain
SST/ELI/ELIB
Bâtiment Carnoy, c.145
Croix du sud 4-5, bte L7.07.04
1348 Louvain-la-Neuve
BELGIUM


	[[alternative HTML version deleted]]