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

Edzer Pebesma edzer.pebesma at uni-muenster.de
Sat May 28 12:15:44 CEST 2016 


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

-- 
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/20160528/3e72af7b/attachment.bin>

More information about the R-sig-Geo mailing list