[geos-devel] rgeos interface to R classes (sp, others)

[Roger Bivand]

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
>> 
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no