[R-sig-Geo] rgeos::gDifference

[Robert J. Hijmans]

You can also do:

library(raster)
x <- sp1 - sp2
y <- sp2 - sp1

(equivalent to the "erase" function)



On Mon, Jan 5, 2015 at 5:43 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> 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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo