[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