[R-sig-Geo] rgeos::gDifference

[Roger Bivand]

On Mon, 5 Jan 2015, Edzer Pebesma wrote:

>
> On 01/05/2015 01:53 PM, Roger Bivand wrote:
>> On Mon, 5 Jan 2015, Edzer Pebesma wrote:
>>
>>> When running the following script:
>>>
>>> library(sp)
>>> gt = GridTopology(cellcentre.offset = c(0,0), cellsize = c(1,1),
>>>     cells.dim = c(3,3))
>>> sgdf = SpatialGridDataFrame(SpatialGrid(gt),
>>>    data.frame(r = round(runif(9, max = 10), 1)))
>>>
>>> gt = GridTopology(cellcentre.offset = c(.25,.25), cellsize = c(1,1),
>>>     cells.dim  = c(1,1))
>>> grd.red = SpatialGrid(gt)
>>>
>>> library(rgeos)
>>> sp1 = as(sgdf, "SpatialPolygonsDataFrame")
>>> sp2 = as(grd.red, "SpatialPolygons")
>>>
>>> plot(sp1, axes = TRUE)
>>> plot(sp2, border = 'red', add = TRUE)
>>>
>>> dif = gDifference(sp2, gUnionCascaded(sp1))
>>> dif = gDifference(sp2, sp1)
>>>
>>>
>>> the last command gives the following error:
>>>
>>> Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_not_poly,
>>> "rgeos_difference") :
>>>  TopologyException: no outgoing dirEdge found at -0.25 -0.25
>>> Calls: gDifference -> RGEOSBinTopoFunc -> .Call
>>> Execution halted
>>>
>>> What does this mean? Why does gUnionCascaded(sp1) not do this?
>>>
>>
>> gUnionCascaded() is deprecated, use gUnaryUnion(). This isn't the
>> problem, the argument order is probably the cause:
>>
>> dif = gDifference(sp1, sp2, byid=c(TRUE, FALSE))
>> length(dif)
>>
>> works, preserving the input geometries, as does:
>>
>> dif = gDifference(gUnaryUnion(sp1), sp2)
>> length(dif)
>>
>> if you want to aggregate.
>>
>> This:
>>
>> dif = gDifference(sp2, gUnionCascaded(sp1))
>>
>> gives a NULL result, as sp2 is within sp1, as does:
>>
>> dif = gDifference(sp2, gUnaryUnion(sp1))
>>
>> but:
>>
>> dif = gDifference(sp2, sp1, byid=c(FALSE, TRUE))
>>
>> preserves the input geometries. If it was clear what you wanted as
>> output from the operation, it would be easier to see what is going on.
>
> Thanks!
>
> I wanted to know if there was any part of sp2 that is not covered by
> sp1, across the entire objects sp2 and sp1. Given the documentation, I
> had expected
>
> gDifference(sp2, sp1)
>
> to give me that, and didn't understand the error message. I'm fine to
> gUnaryUnion() sp1, but don't understand why.

If you don't need the output of the operation, then:

gWithin(sp2, gUnaryUnion(sp1))

or

gCovers(gUnaryUnion(sp1), sp2)

maybe easier. The logic in GEOS as we have it now is that the predicates 
look at the component geometries, it seems, so gUnaryUnion() clarifies 
the query.

Roger

>
> Best regards,
>

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: Roger.Bivand at nhh.no