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

[Roger Bivand]

On Mon, 1 Feb 2010, Martin Davis wrote:

> Can you post the WKT/WKB of your test case?

If this is WKT, then:

> SP2WKT(spol0)
[1] "MULTIPOLYGON(((0 0.1,0 0.2,0.1 0.2,0.1 0.1,0 0.1)),
   ((0.1 0.1,0.1 0.2,0.2 0.2,0.2 0.1,0.1 0.1)),
   ((0 0,0 0.1,0.1 0.1,0.1 0,0 0)))"
> SP2WKT(spolx0)
[1] "MULTIPOLYGON(((0 100,0 200,100 200,100 100,0 100)),
   ((100 100,100 200,200 200,200 100,100 100)),
   ((0 0,0 100,100 100,100 0,0 0)))"
>

I don't think that there is an OGR driver, is there? So this is a text 
representation for the two objects, each with 3 Polygon members, and only 
differing in scale (the coordinates of the second are 1000 times those of 
the first). I see the same difference on GEOS 3.2.0/C API 1.6.0 on x86_64 
RHEL5 and with an MSYS-built GEOS 3.2.0/C API 1.6.0 on Win 32 R.

Hope this is good enough, grateful for any ideas,

Roger


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

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