[geos-devel] rgeos interface to R classes (sp, others)
Martin Davis
mbdavis at refractions.net
Mon Feb 1 13:26:57 EST 2010
Can you post the WKT/WKB of your test case?
Roger Bivand wrote:
> On Mon, 1 Feb 2010, Martin Davis wrote:
>
>> Roger,
>>
>> The technique of using buffer(0) to union polygons is now deprecated
>> in favour of using Unary Union (Geometry.union() - not sure what the
>> exact GEOS signature is). Unary Union is usually faster and more
>> robust than the previous technique. You might want to check this out.
>
> Martin,
>
> Thanks, I've replaced buffer(0) by GEOSUnionCascaded (C API). But the
> artefact remains. For input MULTIPOLYGON:
>
> ---------
> | | |
> | | |
> |--------
> | |
> | |
> -----
>
> I still get the same MultiPolygon out when the bounding box is
> (0,0)(0.2,0.2), but correctly Polygon:
>
> ---------
> | |
> | |
> | -----
> | |
> | |
> -----
>
> when the bounding box is (0,0)(200,200). May this be related to the
> input Polygon objects only having 5 coordinates (4 corners and
> closure)? The generating code is inserting exactly the same doubles
> into the coordinates.
>
> Anyway, I'm sure that Unary Union is a better solution than buffer(0)
> in general, so thanks for that!
>
> Roger
>
>
>>
>> Martin
>>
>> Roger Bivand wrote:
>>> Hi,
>>>
>>> I've felt that I've been making reasonable progress with interfacing
>>> GEOS geometries and methods for the R language (cran.r-project.org),
>>> in a draft contributed package rgeos:
>>>
>>> https://r-forge.r-project.org/projects/rgeos/
>>>
>>> However, I encountered a problem that I do not understand, and would
>>> be very grateful if I could be put back on track. The specific
>>> problem is that an R function uses the C API to dissolve polygon
>>> borders for adjacent polygons sharing a category given as its second
>>> argument.
>>>
>>> I've used the buffer technique from the JTS documentation, and all
>>> was well until I tried to dissolve borders between touching squares
>>> when the coordinate measures were small (square side 0.1). When the
>>> squares are 100 units, all is well, and GEOMTouches is TRUE. But for
>>> 0.1, GEOMTouches is FALSE, and no dissolving takes place.
>>>
>>> This can be reproduced (I'm on RHEL5, x86_64) by installing R,
>>> contributed packages sp and maptools from CRAN, GEOS (3.1.1 or
>>> 3.2.0), and installing rgeos from R-Forge. Start R, say
>>>
>>> library(rgeos) example(unionSpatialPolygonsGEOS)
>>>
>>> and look for undissolved borders in:
>>>
>>> image(grd, axes=TRUE)
>>> plot(spol1, add=TRUE)
>>>
>>> but dissolved in:
>>>
>>> image(grdx, axes=TRUE)
>>> plot(spol1x, add=TRUE)
>>>
>>> I started on an alternative implementation using GEOMTouches in
>>> unionSpatialPolygonsGEOS(..., buffer=FALSE) output to console, where
>>> one sees in the example output:
>>>
>>>> spol1F <- unionSpatialPolygonsGEOS(as(spol, "SpatialPolygons"),
>>> + as.character(spol$xvs), buffer=FALSE)
>>> # 4 squares within (0,0), (0.2,0.2) NE, NW, SW share category, SE
>>> doesn't
>>> npls: 3, nnpls: 3
>>> type[0] Polygon
>>> i: 0, j: 1, touches: 0
>>> i: 0, j: 2, touches: 0
>>> type[1] Polygon
>>> i: 1, j: 2, touches: 0
>>> type[2] Polygon
>>> out of function
>>> npls: 1, nnpls: 1
>>> type[0] Polygon
>>> out of function
>>>> spol1xF <- unionSpatialPolygonsGEOS(as(spolx, "SpatialPolygons"),
>>> + as.character(spolx$xvs), buffer=FALSE)
>>> # 4 squares within (0,0), (200,200), same categories
>>> npls: 3, nnpls: 3
>>> type[0] Polygon
>>> i: 0, j: 1, touches: 1
>>> i: 0, j: 2, touches: 1
>>> type[1] Polygon
>>> i: 1, j: 2, touches: 1
>>> type[2] Polygon
>>> out of function
>>> npls: 1, nnpls: 1
>>> type[0] Polygon
>>> out of function
>>>
>>> The C code is in:
>>>
>>> https://r-forge.r-project.org/plugins/scmsvn/viewcvs.php/pkg/src/rgeos_sp.c?rev=44&root=rgeos&view=log
>>>
>>> around line 466 rgeos_SpatialPolygonsUnion(), calling
>>> rgeos_plspairUnion() - desperate test framework, or
>>> rgeos_plsbufUnion() which I had thought worked, but which clearly
>>> doesn't dissolve small squares.
>>>
>>> Very grateful for any help, this is a show-stopper, and I had hoped
>>> to release rgeos before February (small chance now!)
>>>
>>> Roger
>>>
>>
>>
>
--
Martin Davis
Senior Technical Architect
Refractions Research, Inc.
(250) 383-3022
More information about the geos-devel
mailing list